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Active galactic nuclei, x-ray binaries, pulsars, and gamma-ray bursts are all be- 
lieved to be powered by compact objects surrounded by relativistic plasma flows driv- 
ing phenomena such as accretion, winds, and jets. These flows are often accurately 
modelled by the relativistic magnetohydrodynamics (MHD) approximation. Time- 
dependent numerical MHD simulations have proven to be especially insightful, but 
one regime that remains difficult to simulate is when the energy scales (kinetic, ther- 
mal, magnetic) within the plasma become disparate. We develop a numerical scheme 
that significantly improves the accuracy and robustness of the solution in this regime. 
We use a modified form of the WENO method to construct a finite-volume general 
relativistic hydrodynamics code called WHAM that converges at fifth order. We avoid 
(1) field- by-field decomposition by adaptively reducing down to 2-point stencils near 
discontinuities for a more accurate treatment of shocks, and (2) excessive reduction 
to low order stencils, as in the standard WENO formalism, by maintaining high order 
accuracy in smooth monotonic flows. Our scheme performs the proper surface integral 
of the fluxes, converts cell averaged conserved quantities to point conserved quantities 
before performing the reconstruction step, and correctly averages all source terms. We 
demonstrate that the scheme is robust in strong shocks, very accurate in smooth flows, 
and maintains accuracy even when the energy scales in the flow are highly disparate. 

Key words: accretion, accretion discs - black hole physics - galaxies: jets - hydro- 
dynamics - magentohydrodynamics (MHD) - methods: numerical 



Astrophysical systems containing compact objects are the 
brightest and most efficient engines in the Universe. Enor- 
mous amounts of energy in the form of radiation and out- 
flows are released by accreting blac k hole syste ms at the 
heart of active gala ctic nuclei (AG N, iKrolikll 19991 ). gamma- 
ray bursts (G RBs, Woosley | ll993l ). and black hole X-ray bi- 
naries (XRBs, [L ewin et al.lll995f ). Pulsars and their associ- 
ated plerions are powered by th e relativistic spin-down o f 
highly magnetized neutron stars l|Goldreich fc Julianl fl969). 
These systems involve strong magnetic fields, relativistic 
flows, and strong gravity. A general framework that has 
proven to be accurate in this regime is the ideal ma gne- 
tohydrodynamic approximation (MHD) (Phinncy 1983). 
Significant uncertainty has remained in how black hole 



* E-mail: atchekho@cfa. harvard.edu (AT); 
jmckinney@cfa.harvard.edu (JCM); 
rnarayan@cfa.harvard.edu (RN) 



accretion systems are powered and produce relativistic jets. 
Severa l seminal works have ou tlined plausible key mecha- 
nisms. iBalbus fc Hawlevl!|l99ll ) showed that the magnetoro- 
tational instability (MRI) drives the turbulent transport of 
angular momentum within an accretion disk. The MRI gen- 
erates a dynamo that a mplifies any weak magnet ic field 
in the accreting plasma. iBlandford fc Pavnel l|l982f ) found 
self-similar MHD solutions that describe winds launched 
magnetocentrifugally from thin magnetized accretion disks. 
Blandfor d fc Znaie 

2 (|l977l ) described a model that shows 
how black hole spin energy could be used to power a rel- 
ativistic jet. Despite these important advances, it remains 
unclear how the various mechanisms coexist and how effec- 
tive they are in real astrophysical systems. 

Time-dependent numerical simulations have demon- 
strated that an accreting rotating black hole produces 
a disc wind and a pair of relativistic jets as a natu- 
ral outcome of quasi-s teady magnetized disc accretion. 
|Pe Villiers et all l|2003l ) identified the general morphol- 
ogy of black hole accretion flows to consist of a disk, a 
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corona, a coronal outflow, a nd an evacuated funnel where 
the magnetized jet develops. iMcKinnev fc Gammiel (|2004l ) 
sho wed that the magnet ized funnel jet is well-described by 
the iBlandford fc Znaiek solution. Follow ing up on the pio - 
neering simulation s by lKoide et all (120021 ) and lKoidd (|2003h , 
iKomissarovl (|2005h showed that even without a disk the nat- 
ural o utcome is that black hole spin energy i s extr acted 
via the Blandford fc Znaiekl effect. Iffirose et ail [|2004h and 
IMcKinnev! (|2005l ) showed that the field geometry associ- 
ated with the Blandford-Znajek effect is the most sta- 
ble large-scale field geometry in the presence of a turbu- 
lence disk around a black hole. iKomissarov fc McKinnevI 
<|2007h showed that the Blandford-Znajek effect operates 
even for a = 1 despite concerns that at high blac k hole spin 
the fiel d might get expelled from the horizon. iMcKinnevI 
(2006b) demonstrated that the disc models studied by 
|Pe Villiers et all (|2003h and IMcKinnev fc Gammiel (12004m 
produce a self-consistent magnetized jet with a Lorentz 
factor of 7 ~ fO and an opening angle of ~ 5°, 
demonstrating that the ideal MHD approximation is suffi- 
cient to explain many r elativi stic astrophysical jets. Finally, 
IMcKinnev fc Naravanl ( |2007al lbh showed that the Blandford- 
Payne magnetically-driven wind does not properly describe 
the wind from turbulent disks near black holes and that 
the low density magnetized jet from the black hole is colli- 
mated by the disk corona and wind rather than being self- 
collimated. 

Many of the results described above were obtained 
from simulations that used numerical methods based upon 
so-called 'conservative schemes,' which have been proven 
to be accurate and robust for mo d eling relativistic flows 
jAlov et all Il999l. |2000|; iFontl 120031; iMarti' fc Mulled 120031 ; 
iLeismann et alj|2004 lAlov et al.ll2005l ). The main advan- 
tage of these schemes is that they conserve the integrals of 
motion up to machine precision, and this explicit conserva- 
tion enables them to accurately resolve discontinuities in the 
flow. However, systems with kinetic, thermal, or magnetic 
energy scales differing by orders of magnitude pose serious 
difficulties for conservative numerical schemes. This happens 
in, e.g., highly supersonic jets and winds where the internal 
energy is much smaller than the kinetic energy. Since con- 
servative schemes evolve the values of momenta and total 
energy (the sum of the kinetic, internal, and magnetic ener- 
gies), even a relatively small error in the evolution of these 
quantities can destroy the accuracy of such flows. In the hy- 
drodyna mic case, this is referred to as the h igh Mach number 
problem (iRvu et al.lll993l ; lFeng et al.ll2004) . A similar prob- 
lem occurs in the highly magnetized MHD regime, where the 
equations become stiff when the effective magnetic "Mach" 
number M mag = ^Ja f , where a is the comoving mag- 
netic energy density per unit rest mass energy density. In 
this paper we focus on solving the hydrodynamic high Mach 
number problem, although the ultimate goal is to eventually 
apply our method to the MHD equations to study black hole 
and neutron star systems that may have M mag > 10 3 . 

There are alternatives to evolving the total energy 
equation and some of them provide a reasonably ac- 
curate treatment of high Mach number flows. On e can 
use a ZEUS-type method (|Stone fc Normanl ll992T l that 
solves the internal ener gy evolution e quation with an ar- 
tificia l viscosity (e.g., Hawlevl I2OO0I : lleumenshchev et al.l 
120031 ; |Pe Villiers fc Hawlevl 120031 ). However, the schemes 



by lHawlevI (|2000l ); |Pe Villiers fc Hawlevl (|2003l ) lose en- 
ergy gener ated in reconnectin g curr ent sheets and the 
scheme by lleumenshchev et al.l ((2003), which uses an ar- 
tificial resistivity, still does not strictly conserve energy. 
Bucciantini ct al. I (|2006h employ a method that solves the 
entropy evolution equation and does not allow any shocks 
in the solution. 

Another approach is the 'dual-energy formalism' that 
involves switching from evolving the total energy to evolv- 
ing the internal energy (or entropy) equations of motion 
for a Mach number large r than a th r eshold value M t h- For 
instance, M t h ~ 10 for IRvu et al.l dl993h. Mth ~ 5 for 
iBrvan et all d 19951) . and M th ~ 50 for lFeng et all (|2004h . For 
example, IRvu et al.l (|l993T ) switch from evolving the energy 
to evolving the entropy equation when (1) the Mach num- 
ber is greater than ~ 10 and (2) either the shock strength is 
< 1/3 or the flow is diverging. This criterion can fail to cap- 
ture the flow properly in the presence of weak shocks (which 
could dominate the energy dissipation). Also, once the en- 
tropy equation is used to determine the pressure, the pres- 
sure could be secularly underestimated if a shock is building 
up. 

Finally, there a re Lagrangian-type schemes such as 
bv lTrac fc Peril (2004). who use a moving grid technique for 
treating high Mach number flows (see section [5721 for further 
discussion). However, their scheme requires a sophisticated 
treatment of the advection terms in the equation of motion. 
There also exist sp ecialized schemes that are independent of 
the Mach number (|Marv et al.ll2000[ ). 

We describe a general relativistic hydrodynamic scheme 
that solves the total energy equation yet provides an accu- 
rate and robust solution in the difficult regime of high Mach 
number flows. This is achieved by (1) being finite-volume 
(conserves integrals of motion to machine accuracy), (2) be- 
ing fifth order accurate, (3) maintaining fifth order accuracy 
in monotonic regions (which avoids the excessive reduction 
that occurs in the standard WENO formalism), (4) using 
adaptive stencil reduction near discontinuities, and (5) ob- 
taining fluxes from reconstruction of primitive, rather than 
e.g. conserved, quantities. Our scheme performs the proper 
surface integral of the fluxes and converts cell averaged con- 
served quantities to point conserved quantities before the 
reconstruction step. It also averages the source terms that 
appear on the right hand side of the conservation equa- 
tions. These averaging and de-averaging procedures enable 
our scheme to have an error term that is fifth order in grid 
cell spacing for an arbitrary metric and coordinate system. 
Our scheme use s a modified form of the WENO method (see, 
e.g., [Sh7Jll997h for spati al discretisation and the method of 
lines appr oach (see, e . g.. lLeVeaudll99ll) with Runge-Kutta 
stepping (jPress et al.l Il992l ; IShulll997l ) for time discretisa- 
tion. We eliminate the expensive eigenvector decomposition 
step used by other high order methods by adaptively reduc- 
ing down to 2-point stencils in the vicinity of discontinu- 
ities. This reduction technique gives a similar accuracy to 
eigenvector decomposition without the additional computa- 
tional expense. On the other hand by maintaining high order 
accuracy in smooth monotonic flows, we are able to avoid 
excessive reduction to low order stencils that occurs in the 
standard WENO formalism. 

Our numerical scheme is called WHAM, which stands 
for WENO high accuracy magnetohydrodynamics. In sec- 
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tion [2] we describe the equations that determine the rele- 
vant physics (ideal MHD, single fluid approximation) in the 
form used by the scheme. This is followed by a detailed de- 
scription of the finite- volume numerical scheme in section [3] 
and of the fifth order WENO-type reconstruction procedure 
in section [4] Section [5] presents results of an extensive set 
of tests for code verification. Section [6] describes the lim- 
itations of the numerical scheme, section [7] outlines some 
possible applications, and section [H] summarizes our results 
and concludes. In the Appendix we describe the method for 
reducing the stencil size near discontinuities, preserving the 
fifth order of reconstruction in monotonic regions, and other 
implementation details. 



2 GOVERNING GRMHD EQUATIONS 

While this paper focuses on a new purely hydrodynamical 
numerical scheme, many of the tools we develop can be ap- 
plied to a full MHD scheme. Therefore we keep our dis- 
cussion as general as possible. Here we outline the general 
relativistic magnetohydrodynamic (GRMHD) equations of 
motion in the form of conservation laws. Throughout, we 
use the Einstein summation convention, where Greek let- 
ters run from to 3 and Latin letters from 1 to 3. The units 
are such that G = c = 1. 

The continuity equation that describes baryon number 
conservation is given by 

if**)* = 0, (la) 

where p is mass density in the fluid frame, u M is the 4- velocity 
of the fluid and the '; p' subscript denotes the covariant 
derivative with respect to . The energy-momentum con- 
servation equation is given by 

= 0, (lb) 

where T Ml/ are the components of the energy-momentum 
tensor. Finally, the source-free part of Maxwell's equations 
describes the evolution of the fields, 

= 0, (lc) 

where *F M " = 1 / 2 t^ VKX F K x ar e the components o f the dual 
of the electromagnetic tensor (|Misner et al.lfl973h . 

The energy-momentum tensor can be written as a sum 
of matter and electromagnetic contributions, T^ v — T^" + 
Te'm , with 

T^ = {p + u g +p g )uV+p g g^, (2) 
= F^F\ - 7 4 g^F KX F nX , (3) 

where u g and p g are the internal energy and pressure in 
the fluid frame, and F^ v is the electromagnetic field ten- 
sor. Introducing the comoving magnetic field 4-vector, V = 
u lL *F y,v , in ideal MHD it can be shown that 

T\ = (p + u g + Pg + b 2 )u"u v + (p g + 7 2 b 2 ) b\ - b"b v (4) 

(|Gammie et al.ll2003h . 

For numerical purposes we write down equations (|la[) 
- (|lc[) in a coo rdinate basis. The contin uity equation (|la[) 
takes the form (|Landau fc Lifshitz|[l975h . 

ftCV-PPU*) + dii^/^pu) = 0, (5a) 
where g = Det g^ v is the metric determinant and the 't' 



index denotes the time component. Similarly, the energy- 
momentum equation ()lb|) takes the form 

at(V=?T*„) + diW=gi* v ) = V=gT\r\ K , (5b) 

where r VK are the connection coefficients. Finally, with 
the introduction of 3- vectors of magnetic and electric fields, 
B l = *F U and E l = F"/\/— ff; equation |lc]> is equivalent to 
the induction equation, 

d t (V^B i ) = -d j [V^(e» k E k )] 

= -d j [V=9(B i v j -B j v% (5c) 

where if — u l /u is the 3- velocity, plus the no-monopoles 
constraint, 

9 i (y/=gB i )=0. (5d) 

In deriving the last equality in equation (|5c|l we have used 
the ideal MHD approximation Ei — —etjkv-' B = — (vxB)j. 

Equations (|5a|| - ()5d|) determine the evolution of a 
plasma system for given initial and boundary conditions. 
We derive numerical discretisations of these equations in the 
next section. 



3 NUMERICAL SCHEME 

We develop a numerical scheme based upon a solution to a 
general set of conservation laws. 

3.1 Vector form of conservation laws 

The equations of motion (|5a[l - (|5c[l can be written in the 
form of a vector conservation law: 

a t U(P) + 9 I F 1 (P) = S(P), (6) 

where the vector of conserved quantities is 

U(P) = y^(p M t ,T t t +p M i ,T t J ,B fc ), (7) 

the vector of fluxes in the ith direction is 

F'(P) = ^-g{pu\T\ + pu\ T i i ,B t v k - B V), (8) 

and the vector of source terms is 

s(p) = v=i(o, T" A rV, o, o, o), (9) 

where, according to the convention, k run from 1 to 3 
and A, p, v run from to 3 so the vectors U, F ! , and S have 
8 components each. We have analytically removed the mass 
energy density from the conserved quantity that corresponds 
to the total energy, the second component of U, so that 
the method remains accurate in the limit of nonrelativistic 
flows (see Appendix lA9|l . As a set of primitive quantities we 
choose 

P = (p,ug,u j , B k ), (10) 

where the spatial components of the relative 4-velocity 
v? = u 3 — jr/ 3 , in which 7 = —u^r}^ is the Lorentz factor 
of the flow as measured in the normal observer frame, and 
?7 M = (—a, 0, 0, 0) is the 4-velocity of a zero angular mo- 
mentum observer in the coordinat e basis for axisym metric 
space-times; here a = (-g u )~ 1/2 l|Noble et all 12006). One 
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can show that u 1 = and 7 = yl + Wui. Using the relative 
4- velocity (as compared to the usual 4- velocity) has 2 advan- 
tages: (1) its interpolation always leads to a physical state 
and (2) even in rotating (e.g., Kerr) space-times its spatial 
components determine the unique value of w* = J /a. 

Note that system © does not contain the no-monopoles 
constraint (|5dl) . This constraint is considered in future work, 
see section [6] 



3.2 Motivation for a consistent finite volume 
scheme 

Eulerian conservative numerical schemes keep track of the 
values of conserved quantities for each grid cell. In this sec- 
tion we show the importance of accounting for the difference 
between the cell averaged and cell centred values of these 
quantities for maintaining high accuracy in studies of highly 
supersonic hydrodynamic flows. In such flows small relative 
errors in the total energy can lead to large relative errors in 
the internal energy. This is known as the high Mach num- 
ber p roblem <|Rvu et all 1 19931 : iTrac fc Penll2004l : iFeng et all 
2004). We show that even though the errors are second order 
in grid cell spacing, their magnitude is proportional to the 
square of the Mach number and so these errors can destroy 
the accuracy of numerical models. 

Consider a numerical method based upon piecewise lin- 
ear interpolations of primitive quantities, and consider a flow 
in which p oc const, v oc x, and the internal energy is neg- 
ligible relative to the kinetic energy, u g <C pv 2 . Is a linear 
interpolation of the conserved energy then sufficient? That 
is, can we assume that there is a negligible difference be- 
tween E(x) = pv 2 /2 + u g oc x 2 at the centre x = Xi of a 
grid cell, and its average value (E) over that grid cell? The 
answer is we cannot since for E{x) defined above, 



{E}-E = % E"( Xl )Ax 2 = 7: 



24 PV 



> U„ 



(11) 



where Ax is the grid cell spacing. Neglecting the difference 
between E and (E) will manifest as a significant error in the 
smaller quantity u g . 

In particular, consider a simple nonrelativistic one- 
dimensional, high Mach number hydrodynamic flow (no 
gravity or any other external forces). Suppose initially the 
flow has a uniform density and pressure, and its velocity is a 
linear function of position. We call this a Hubble-type flow. 
The primitive flow variables (density, velocity, and specific 
internal energy) at time t are given by: 



Po 



P(M) = TT^P' 
«(*,*)- v '° x 

Ug{X,t) 



(12) 



Uo 



(i + v> ty> 

where v' = dv/dx(x,t — 0). We assume the equation of 
state p g = (r — l)w B , where F is the adiabatic index of the 
gas, and we choose the flow parameters po, uo, and v'o such 
that the flow is highly supersonic in the grid cell centres, 

Ug < PV 2 . 

As we show in Appendix lAll if one discretises the above 
equations in space and time to any order but neglects the 
difference between E and (E), then after one time step one 




i+I/2,} 



Figure 1. Location of quantities within the grid cell Ay. Thick 
solid lines show the grid of cell interfaces and dashed lines show 
the grid of cell centres. The cell centre is indicated by the large 
dot where the primitive P{j, point conserved Uij, and average 
conserved (U)ij quantities are located. The fluxes are located 
at the centres of grid cell interfaces shown by smaller dots. See 
sections 13.31 - 13.41 



makes a relative error in internal energy, 
^1 ~_Af 2 - St 

Uo 



(13) 



where M^ in ~ pov'o Ax 2 /uq is the square of the minimal 
value of the Mach number on the grid, Aa; is the grid cell 
size, and 8t — v' At is the dimensionless time step for the 
problem for a computational box with outer edge x = 1. 

At every time step the internal energy uniformly de- 
creases by more than it should and after several time steps 
the error will be quite large, and u g may even become nega- 
tive. Even though this error in the internal energy is second 
order in space, the coefficient of this error is proportional 
to the Mach number of the flow squared. Therefore, for a 
large Mach number, any scheme that does not distinguish 
between cell averaged and point conserved quantities has to 
use a resolution proportional to the Mach number of the flow 
in order to correctly capture the evolution of the internal en- 
ergy over the relevant timescale of 5t ~ 1. For instance, for 
Mmi n = 100 such a scheme requires roughly a lOOx resolu- 
tion increase. See section 15.11 for a numerical verification of 
these statements. 



3.3 Numerical grid 

The numerical scheme is built upon a uniform grid in a 
coordinate basis, where for simplicity we consider here the 
two-dimensional case (the three-dimensional case is a trivial 
extension - also see section [BJ. The grid consists of cells 



that we define as Aij — (x]_ 



1/2! • 



2) x 



-l/2i Jj j+1/2J^ 



where {x\ , x 2 ,t n ) = (iAx 1 , jAx 2 , nAt), Ax 1 and Ax 2 are 
the cell spacings (see figureflj, and At is th e time step. For a 
discus sion of numerical grid generation, see lThompson et al.l 
1 19851 ). 
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An arbitrary smooth transformation can be used to map 
this computational grid into physical space as a fixed adap- 
tive mesh that focuses on the regions of interest (for more 
detail, see lGammie et al"]|2003l ). The transformation of any 
tensorial quantity requires derivatives of the coordinate map 
that are either computed analytically or using a simple nu- 
merical scheme that results in an ac curacy of roughly the 
machine accuracy to the 2/3 power (|Press et al ] l!992fl . A 
similar method is used to compute the connection coeffi- 
cients. Therefore the metric identities are of similar accu- 
racy. Both the discretisation of equations of motion and the 
reconstruction of functions are independent of the physical 
coordinates, and this makes the algorithm simple and gen- 
eral. 



3.4 Finite volume discretisation of the equations 

In this section we derive a finite volume discretisation of the 
conservation equations (|6]) for a two-dimensional problem. 
This discretisation guarantees that the integrals of motion 
are conserved up to machine precision error. Integrate any 
component of the vector equation ([6]) over a grid cell V = 
Aij to obtain 



dx 



JJ^ dxldx2 = JJ Sdxl 



dx 



(14) 



Taking the integrals and dividing by grid cell volume 
Ax 1 Aa; 2 , we obtain a syst em of semi-d iscrete equa- 
tions (method of lines, see e.g. iLeVcauc 199j]): 



d(JJ)n , {F 1 ) i+1 / 2 ,j - (J ?1 ) i -i/ 2 , J 



dt 



+ 



Ax 1 



(^ 2 k J + l /2 - {F 2 )j 



ij-l/2 



= {S)a, (15) 



where the angle brackets denote spatial averaging over the 
volume/surface of the grid cell for quantities located in- 
sidc/at the surface of the grid cell (Figure [TJ. So far this 
equation is exact. 

We use R unge-Kutta 4th order accurate time step- 
ping (RK-4, see lPress et aHll992T ) for discretising the system 
of ordinary differential equations (|15[) in time. Each Runge- 
Kutta trial step is a first order numerical discretisation of 
equation (|15|) given by: 



11,3 



At 



+ 



(F 2 ) 



Aa; 1 
~ {F 2 } 



Ax 2 



= {S)a. (16) 



The RK-4 algorithm combines a series of four first-order 
trial time steps, which we refer to as substeps, with different 
At's to obtain a 4th order accurate analog of (|16[) . We note 
that formally we need to use RK-5 time stepping in order 
for our scheme to converge at fifth order. However, for most 
problems RK-4 produces time stepping errors that are much 
smaller than the truncation errors due to spatial reconstruc- 
tion, so RK-4 is practically suffic ient for attaining fifth order 
convergence (c.f. section [Ol and IZhang fc MacFadvenll2006l . 



hereafter IZMOd : however, see section[5TT} ■ We find that RK-4 
alway s gives mo re accurate results than the RK-3 scheme 
from[Shu| l| 19971 ). for both smooth and discontinuous flows. 
We solve each Runge-Kutta substep given by equation (|16|l 
using the Godunov technique. 



3.5 High-order Godunov schemes 

Godunov schemes perform the evolution in time of a set of 
grid cell averages of conserved quantities, {(U)^}, by con- 
sidering binary interactions between adjacent cells. These in- 
teractions are usually approximated to be one-dimensional. 
In the simplest case, the first order Godunov scheme, the 
distribution of conserved quantity inside each grid cell is as- 
sumed to be a constant. The one-dimensional interaction be- 
tween two such distributions is the classical Riemann prob- 
lem: the d ecay of a d iscontinuity between two constant dis- 
tributions l|Toro|ll997M . For a given discontinuity, an exact or 
approximate Riemann solution is obtained for the flux used 
to update the conserved quantity and obtain {(U}™ +1 } at 
t — t n +\. These first order schemes are robust but inaccurate 
for smooth flows or flows with discontinuities not modelled 
by the Riemann solver. 

Higher-order Godunov schemes attempt to obtain 
higher accuracy by using higher order methods to com- 
pute both the interface states and the temporal updates. 
High-order methods are desirable since typically their com- 
putational cost is much less compared to increasing the 
numerical resolution for first order schemes. Higher-order 
Godunov numerical schemes generate a discontinuity at 
each interface as a result of a high-order reconstruction 
within adjacent grid cells (see section [4|. The left /right in- 
terface states are fed to the Riemann solver that is used 
to obtain the flux at the interface. However, such Rie- 
mann solutions actually assume the distributions are con- 
stant within each cell and that the characteristics of the 
Riemann fan are linear ( see, e.g., Colella fc Woodward! 
1 19841 : iMignone et ai1l2005l : IZhang fc MacFadvenll2006h . For 
example, one of the most advanced reconstructi on schemes 
is the PPM scheme (IColella fc Woodward1ll984r i where the 
left/right states (fed to a standard Reimann solver) are ob- 
tained by an average over a 'domain-of-dependence.' This 
averaging partially accounts for the left/right distributions 
being parabolic, but it still assumes the characteristics are 
linear. 

Since high-order schemes actually generate non- 
constant distributions in each cell, the Riemann fan is non- 
linear and so standard Riemann solvers give an incorrect 
flux at the interface. This inconsistency leads to incorrect 
eigenwave amplitudes and potentially to spurious oscilla- 
tions. In order to compute a Riemann solution consistently 
with the high-order scheme, one should instead solve the 
Generalize d Riemann Problem (GRP, Section 13.4.1 from 
lTorolll997h that involves non-constant background distribu- 
tions and so curved charac teristics for the Riemann fan (see 
also iToro fc Titaxevll2u06h . However, using standard Rie- 
mann solvers may be sufficiently accurate for most prob- 
lems. For simplicity, like most other schemes, WHAM uses 
such a standard Riemann solver. 
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3.6 Algorithm outline 

We now describe our implementation of a consistent finite 
volume scheme, where by consistent we mean the numeri- 
cal method takes into account the differences between point 
and average conserved quantities, fluxes, and source terms. 
Performing the conversion between cell averaged and cell 
centred values of conserved quantities is crucial for any test 
problem with disparate energy scales, e.g. high Mach num- 
ber flows. In sections 13.21 15.11 15.21 15.31 and 15.101 we give 
examples of highly supersonic flows that illustrate the im- 
portance of accounting for this difference. 

The general procedure of our method is outlined below 
and can be used with any definition of the cell centre to 
interface reconstruction, the averaging and de-averaging of 
conserved quantities, the averaging of fluxes, and the averag- 
ing of source terms. Our particular reconstruction methods 
are described in section [4] 

We subdivide the grid into three domains, denoted by 
Roman numerals in figure [2] In the standard boundary do- 
main I the primitive quantities are specified by boundary 
conditions and in the standard computational domain III 
the conserved quantities are evolved. We also introduce an 
additional intermediate domain II where both the conserved 
quantities are evolved and the values of primitive quanti- 
ties are specified according to the boundary conditions. This 
technique allows our boundary routines to work only with 
the primitive quantities in both the original boundary do- 
main I and the intermediate boundary domain II. Without 
this additional domain our scheme would require applying 
the boundary conditions on the conserved quantities as well. 
Thus, the introduction of domain II simplifies our boundary 
condition routines, including thos e that involve parallel com- 
putation. As in the HARM code l|Gammie et al.ll2003h . our 
scheme is parallelized using domain decomposition within 
the Message Passing Interface (MPI) . For reasonably chosen 
resolutions per CPU (e.g. 32 2 to 64 2 ) we achieve no less than 
70% parallel efficiency for 256 CPUs for 2D problems on a 
cluster with dual 2.0 GHz Opteron processors connected by 
Gigabit ethernet. 

The initial conditions are set in domain III (active 
grid cells), where we define primitive quantities at cell cen- 
tres, {Pij}. The boundary conditions for primitive quanti- 
ties are set in domains I & II (boundary grid cells). In do- 
mains I - HI, we then convert the point primitive quantities 
to point conserved quantities {U°j}. In domains II & III, 
the conserved quantities are either numerically or analyti- 
cally averaged over grid cells to obtain {(U)^}. 

The cell-averaged conserved quantities are evolved in 
time in domains II & III by a sequence of Runge-Kutta sub- 
steps given by (|16fl • For each Runge-Kutta substep the cell- 
centered primitive quantities {P™ +1 } and cell-averaged con- 
served quantities {(U)™^ 1 } at the new time t n +i are found 
from previously computed cell-centered primitive quantities 
{P™} and cell-averaged conserved quantities {{U)*j} (de- 
fined to be at time t n ) by the below procedure: 

Step i. Using the cell centre to interface reconstruction 
on cell centred values of primitive quantities {P^}, ob- 
tain their values on both sides of every cell interface in do- 
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face (i + 1 / 2 , j); 



i+l/2-0,j 



and P" 



for the inter- 



Figure 2. Computational grid for an illustrative 4x3 resolution. 
Grid cells are shown as squares with cell centres being the centres 
of these squares. Upon initialization we set the values of primi- 
tive quantities in the white grid squares (domain III, active grid 
cells) according to the initial conditions and in the dark and light 
shaded squares (domains I and II, the boundary cells) according 
to the boundary conditions; additionally, in domains II and III 
we compute the average values of conserved quantities either an- 
alytically or numerically. For domains II and III, at every time 
step we interpolate the primitive quantities from cell centres to 
cell interfaces to obtain the inter-cell fluxes, average them over 
the cell faces, and update the average values of conserved quanti- 
ties in these domains. Finally, in domain III, we obtain the values 
of primitive quantities by converting the average values of con- 
served quantities to point values. In domains I & II the primitive 
quantities are set using the boundary conditions. The size of the 
white active grid cells region can be arbitrarily large (only limited 
by memory and speed considerations), whereas the width of the 
boundary cell layer surrounding it (shaded regions) only depends 
on the order of the scheme used. 



Step i i. Using the app roximate HLL Riemann 
solver (iGammie et alj [2003), obtain the flux at the 
centre of each interface in domains II & III, e.g., 

Tjil (pn pn \. 

r i+l/2,j( tr i+l/2-0,j i r i+l/2+0,jJi 

Step iii. Spatially average the flux through every interface 
in domains II & HI, using the centre to average reconstruc- 
tion, to obtain, e.g., (F 1 ) i+1/2 , J ; 

Step iv. Find the point values of source terms in do- 
mains I - HI, Sij, according to equation §9$; 

Step v. Average the point values of source terms over grid 
cells, using the centre to average reconstruction, to obtain 
the cell-averaged values of source terms in domains II & III, 

Step vi. Compute the cell-averaged conserved quantities 
at the new time in domains II & III, {{U)™ +1 }, by equa- 
tion (fTo) ; 

Step vii. De-average the conserved quantities, using the 
average to centre reconstruction, to obtain the point val- 
ues {U™^~ }, located at grid cell centres, in domain III; 

Step viii. Obtain the corresponding set of primi- 
tive quantities jPJV 1 " 1 1 using a primitive variable 
solver iMignone fc McKinnevI i|2007h ; 

Step ix. Apply the boundary conditions on the primitive 
quantities to get their values at grid cell centres in do- 
mains I & II; 

We refer to the above algorithm as the WHAM scheme. 
For the purpose of demonstrating the importance of dis- 
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tinguishing between average and point values, we also con- 
sider an inconsistent scheme WENO-IFV that is the same 
as WHAM except we disable the averaging and de-averaging 
procedures (in the initial conditions and in |Step iii| |Step v| 
and |Step vii[ ). We also consider yet another scheme, denoted 
as the WHAM-IS scheme, for which we disable only the 
source integration by skipping |Step v| This enables us to 
determine the effect of not averaging the source termsQ 

A feature of our finite volume approach is that con- 
served quantities are conserved to machine precision for all 
time. An alternative is to use the finite difference approach 
that evolves the cell centred point values of conserved quan- 
tities, so that the c onserved q uantities are then conserved 
to truncation error (|Shulll997r i. While the finite difference 
approach can be made to account for the difference between 
averaged and point values of conserved quantities, we find 
that it is less robust than the finite volume method we have 
described and less desirable since it does not conserve the 
proper quantities. 



4 RECONSTRUCTION 

In this section we explain the method WHAM uses to per- 
form the cell centre to interface reconstruction in |Step i] the 
averaging in |Step iii| and |Step v| and the de-averaging in 
|Stcp vii"] 

4.1 Quantities to reconstruct 

High-order schemes are built upon reconstructing some set 
of quantities from the cell centre to the cell interface. Com- 
monly either the primitive or conserved quantities are chosen 
to be interpolated to obtain the interface state. We perform 
the cell centre to interface reconstruction on density and in- 
ternal energy, and for the reconstruction of velocity we use 
a pro cedure that is more robust in relativistic shocks (see, 
e.g., IZM06h . We separately reconstruct the components of 
the relative 3-velocity, v r = m*/7j and 7. We then multiply 
the interpolated value of v l by that of 7 to get the recon- 
structed value of the relative 4-velocity. 

We perform the averaging and de-averaging reconstruc- 
tions component-by-component on the conserved quantities 
and fluxes. To make the evolution of our scheme failsafe, 
especially since an arbitrary interpolation of the conserved 
quantities and fluxes can lead to an unphysical state, we im- 
plement some features that are more likely to give a physical 
state (see Appendix I A8|l . 

4.2 Multi-dimensional reconstruction 

To perform reconstructions in more than one dimension, 
we use the dimension by dimension approach. It uses one- 
dimensional reconstructions, described in the rest of this 
section, as building blocks for a multidimensional recon- 
struction. This method was found to be faster than a 
genuine multi dimensional reco nstruction and have a sim- 
ilar accuracy (|Shi et al ] l2002h . The generalization of the 

1 'IFV stands for Inconsistent Finite Volume, and 'IS' stands for 
Inconsistent Source. 



cell centre to interface reconstruction to multi dimensions 
is straightforward: we perform one-dimensional cell centre 
to interface reconstructions independently for each dimen- 
sion. In several dimensions the average to centre and cen- 
tre to average reconstructions are performed by applying 
one-dimens ional recon structions sequentially dimension by 
dimension We symmetrize this procedure by 

averaging over all possible permutations of the sequences 
for which the one-dimensiona l reconstruction can be ap- 
plied (see also I Alov et al.lll999T l. This allows us to maintain 
the symmetry of the reconstruction and is our default choice 
for this paper. In particular, we maintain exact symmetry 
for two-dimensional problems^ For multidimensional prob- 
lems that do not require the preservation of symmetry, com- 
putational efficiency can be improved by using Strang-type 
splitting. Without loss of accuracy, we could do this for the 
two-dimensional problems described in sections 15.211 15.221 
and !5.23l However, we use the same settings throughout the 
paper for the sake of showing the ability to use WHAM with 
one single set of parameters for a wide range of problems. 

4.3 One-dimensional reconstruction 

Consider a one-dimensional grid along the x axis consisting 
of grid cells A; = {x i _ 1 / 2 ,x i+1 /2) with cell centres located at 
points Xi — xq + hi, where h is the grid cell size. Godunov- 
type schemes require the conversion of quantities from one 
type of discretised representation on the grid to another: 
e.g., from cell-centred values of density pi to cell-interface 
values Pi+x/2- 

Since in this example the particular dependence of p(x) 
is unknown (only known through its discretisation at cell 
centres pi), in order to get the discretisation of density at 
cell interfaces, we first reconstruct a smooth density profile 
inside each of the grid cells. For this, we use the values of pi 
in several grid cells around the grid cell in which the recon- 
struction is being performed; this set of grid cells is called the 
stencil. Then, we combine the reconstructed profiles inside 
each of the grid cells to obtain the global reconstruction p(x). 
Now we can easily obtain the cell-interface (or any other) 
representation of density by evaluating p(x) at the locations 
of the interfaces (or any other locations). By construction, 
the global reconstruction is smooth inside each grid cell and 
is in general discontinuous at cell interfaces. For this rea- 
son, we adopt a convention by which Pi+1/2 = p( 2; i+i/2-o) 
and Pi-1/2 = p( x i-i/2+o)y he. the one-half in the indices is 
replaced with a value infinitesimally smaller than one-half. 

Within piecewise smooth functions the reconstruction 
is not unique, and the goal is to find such a reconstruction 
out of many possible ones that would give high order of 
convergence in smooth regions and lead to minimal spurious 
oscillations near discontinuities. 

2 The internal registers are often higher precision than the pre- 
cision of variables (e.g., 80-bit for Intel 32-bit machines vs. 64-bit 
for double precision variables) so that one must use a compiler 
option (e.g., -pc64 for Intel compilers) to disable the use of such 
internal precision. These compiler options are only used on those 
parts of the code that involve operations with multiple points at 
the same time. In the case of Intel compilers, one must use the 
-mp compiler option as well in order to restrict the optimizations 
to maintain the order of operations. 



© 2006 RAS, MNRAS OOO.HTMl 



8 A. Tchekhovskoy, J. C. McKinney, & R. Narayan 



Along with the cell centre to interface reconstruction 
considered above, this paper makes use of other reconstruc- 
tion types that can be formulated similarly. Section [4.61 be- 
low gives more detail. 

In the regions where the flow is smooth the reconstruc- 
tion should be of high order to ensure fast convergence to 
the actual solution when the grid is refined. However, next 
to a kink or a shock, the scheme should avoid using stencils 
that contain this non-smooth feature and should instead use 
stencils that reside entirely within the smooth regions of the 
flow: failure to do so would lead to spurious oscillations. 

4.4 Overview of existing algorithms 

In this section we review some existing methods for recon- 
struction. The goal of such methods is to provide higher- 
order accuracy than a piecewise constant reconstruction, 
while avoiding spurious oscillations in the presence of dis- 
continuities. 

Common existing approaches are the following: 
MINMOD, Monotonized Central (MC), van Leer, 
Piecewise Parabolic Method (PPM), Essentially Non- 
Oscillatory (ENO) and Weighted Essentially Non-Oscilla- 
tory (WENO) schemes. Throughout the paper, the order 
of the schemes/reconstructions refers to the order of their 
truncation error. 

In general, MINMOD, MC, and van Leer schemes are 
second order in smooth regions of the flow. Parts of the re- 
construction algorithm used in Piecewise Parabolic Method 
(PPM) are 4th order. PPM is a computationally fast re- 
construction method and is able to reso l ve contact discon- 
tinui t ies we ll (Colel la fc Woodward 1984; G ardiner fc Stond 
120051 : iMimone et alJl2005l ). The standard PPM scheme is 
dimensionally split, resulting in second order accuracy even 
for smooth monotonic flows. Also, when used for the recon- 
struction of primitive quantities it is in general second order 
in smooth monotonic regions because it does not account for 
the differe nce between point and average v alues of conserved 
quantities ijBlondin fc Lufkin|[l993l . IZM06h . 

The schemes in the above paragraph provide a total 
variation diminishing (TVD) reconstruction, i.e. they do not 
increase a measure of t he number an d magnitude of func- 
tion extrema with time (|Hartenlll983h . This produces non- 
oscillatory and robust solutions. However, being TVD also 
means these reconstructions satisfy the maximum principle 
that the reconstructed value at an interface lies in between 
the values at the two neighboring cell centres. Therefore, 
near extrema, where the derivative changes sign, the so- 
called clipping phenomenon may occur: extrema are flat- 
tened thereby increasing the truncation error of the scheme 
near smooth extrema and lowering the order of the scheme 
there to first order ()Harten et al.l 1 19871 : Ijiang fc Tadmorl 
1 19981 ). 

Essentially Non-Oscillatory schemes provide a uni- 
formly high order accurate reconstruction, do not suffer from 
the clipping phenomenon, provide a total variation bounded 
reconstruction (TVB ), and give robust solutions for flows 
with discontinuities l)Harten et al.lfl987h . They choose the 
smoothest stencil out of a set of possible stencils of a fixed 
length. The weakness of this approach is that one needs 
to choose one stencil even if there are several equivalently 
smooth ones, e.g., when the function is a constant. There- 



fore, unless certain measures are taken (see, e.g.. IShu1ll997l ). 
numerical noise, which is always present, may start influenc- 
ing the stencil choice procedure. This is referred to as free 
stencil adaptation and may lead to the numerical noise being 
amplified. Further, these schemes require a computationally 
intensive eigenvector decomposition to minimize the Gibbs 
phenomenon near shocks. 

Convex Essen tially Non-Oscillatory schemes (CENO, 
iLiu fc Oshe"r] [l998) are similar to ENO schemes in that they 
choose one out of the possible stencils for performing the 
reconstruction, but include a mechanism for reducing the 
stencil size in discontinuities. The latter feature allows these 
algorithms to avoid the use of expensive eigenvector decom- 
position, making it easier to ap ply them to gene r al rela tivis- 
tic problems, as wa s done by lAnderson et al.l l)2006h and 
iMizuno et all (2006), who achieve third order convergence 
in smooth flows. However, the CENO method can give non- 
smooth results since the stencil choice may change discretely 
between adjacent grid cells. 

Weighted Essentially Non-Oscillatory schemes 
(WENO) are designed to provide smooth numerical flux and 
do no t suffer from the free stencil adaptation problem (|Shul 
Il997h . Instead of using only the smoothest stencil, they 
take a linear combination with coefficients, called weights, 
of all possible stencils of a given size. The weights in the 
linear combination add up to unity and are distributed in 
such a way that stencils that contain discontinuities get 
extremely small weight. Further, the weights are designed 
in such a way that when the function is smooth in all 
stencils, the weights become close to the optimal ones so 
that the resulting linear combination of the stencils gives 
a higher order approximation - the same order as the one 
that the larger, combined, stencil would give. WENO-based 
schemes are advantageous because they are able to both 
capture shocks and accurately resolve complex smooth 
flow structure. Several groups have had success in devel- 
oping WENO-basc d schemes in app l ication to relativistic 
astro physics (| Rahman fc Moore 2005; Zhang & Mac Fadvenl 
I2006T ) and cosmology dFeng et all 120041 : lOiu et al l 120061 ). 
The high accuracy of WENO schemes enables them to 
be us ed in studies of high Mach number astrophysical 
flows (|Haet alj|2003 ). 



4.5 Weighted essentially non-oscillatory schemes 

Because of the advantages of the WENO method, we have 
implemented this scheme in our code. Consider the one- 
dimensional grid along the x axis that was defined in sec- 
tion [473] Its cell centres are located at points Xi = xo + hi, 
where h is the grid cell size. The right interface of the ith 
cell is located at %+i/2 = Xo + h(i + 1/2 — 0). 

We first consider the cell centre to right interface recon- 
struction. Given the values of a piecewise smooth function 
v(x) at the cell centres, w, = v(xi), we aim to reconstruct the 
function value Vi+1/2 at the interface x = x i+1 / 2 - Consider 
k candidate stencils, each of length k: 



S) = b 



Ci-r + fc-l 



}, r = 0, 



(17) 



Each of these stencils produces its own reconstruction of the 
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function value v i+1 / 2 |Shulll997T l: 

fc-i 

v i+i/2 = E Cr i v i-r+j , r = Q,...,k — l. (18) 
i=o 

with truncation error of 0(h k ) (see AppendixlA"2lfor numer- 
ical values of the coefficients dj)- 

Within the WENO scheme, the reconstructed value of 
v(x i+1 /2) is written as a linear combination of the recon- 
structions due to each of the stencils s r : 



(s r ) 
U i+l/2' 



(19) 



the weights u r are all nonnegative and sum up to unity. The 
heart of the scheme is the proper choice of these weights. 

In the regions where v(x) is smooth, we define a 
WENO-ra scheme as one that delivers a reconstruction of 
order n = 2k — 1, the same order as that of the reconstruc- 

*U *.(<): 



„( s ) 



tion f)t'i/2 due to the combined stencil 5* 



v (S) 

U i+l/2 



i+fe-1 

E * 

=i-k+l 



(20) 



with truncation error 0(h n ). Usually, such constants d r 
can be found that for ui r = d r the reconstructions 
given by equations (|19[1 and (|2Up are exactly equiva- 



lent, i.e. v, 



'1/2 



(S) 

«.,,„. We refer to the constants d r as 



the optimal weights; their values are given in Appendix I A2I 
However, the exact identity is a too strict restriction on 
the weights. In fact, the weights u) T can deviate from the 
optimal weights d r in smooth regions and still produce a 
truncation error of the same order: by (|19p . requiring that 



u r = dr + G(h ), r = 0, 



(21) 



guarante es the sam e order of truncation error (cf. formula 
(2.57) in lShulll997l ). Equation (2lJ is the requirement on 
the weights in regions where v(x) is smooth. When there 
are shocks, the stencils that contain discontinuities should 
be given extremely small weight. 



4.6 Types of reconstruction 

We have just outlined the cell centre to right interface re- 
construction. It inputs a set of point values at the centres of 
grid cells, {vi}, and outputs a set of values at the grid cell 
interfaces, {v i+l / 2 }- 

Analogously, we introduce the following types of recon- 
struction: centre to left interface, {vi} — > {1^-1/2}, centre 
to average, {vi} — > {{«)*}, and average to centre, {(«)»} — > 
{vi}. Here {v)i = hT 1 J* A . v(x) dx, A* = (aJi_i/2, Xi+1/2), is 
the average of v{x) over the ith grid cell. Apart from the 
trivial change in the input and/or output not at ion formu- 
lae (|18p - 1)211) hold for all these types of reconstruction as 
well. Appendix l A2I gives the values of the coefficients dj and 
d r for each reconstruction type used by the scheme. 



3 For instance, for the centre to average reconstruction one uses 
(v)i for the symbol of the output quantity instead of v i+1 /2- 



4.7 WENO prescription for weights 

In the WENO scheme l| Jiang fc Shulll996l : Ishulli997f l. the 
weights are computed based on smoothness indicators, (5 r . 
These indicators give a measure of the variation of p r (x) - 
the reconstruction polynomial due to the stencil s r (i) - 
within the grid cell Ai. The smoothness indicators are de- 
fined as: 



1 r x i+i/2 



E 



d n p r (x) 

dx n 



h' 



dx. 



(22) 



where the factors of h are included to remove any grid cell 
size dependence of /3 r . According to the definition (122[) . a 
smoothness indicator is the average over the cell of interest 
of the sum of the squares of all derivatives of the interpo- 
lating polynomial. Therefore, a smoothness indicator mea- 
sures how deviant from a constant the reconstruction is. 
The smaller the smoothness indicator, the flatter the recon- 
structed profile, i.e. the smaller the change of the interpolat- 
ing polynomial p r {x) over the cell of interest. The smooth- 
ness indicators defined by (|22p are specified for WENO-5 by 
equ ation (IA18P in appe ndix I A4I 

IJiang fc Shul (1996) define the unnormalized weights in 
the following way: 

*=<IW < 23 > 

which, when normalized, become: 

LO r = u5 r /f2, (24) 

where = X^r=o ^ r ^ s * ne sum °^ unnormalized weights and 
e is a small positive parameter. 

The e parameter in (I23p is introduced to avoid divi- 

cilrtn r rv i\ / ^ \ T a 1 n j^ti 1 linn sA r iirAi-l^Ai^t' Trn r-fT -pT-^m 1 H ^ 



10 -6 Oiang & Shul 1996I: 


Shul Il997l 


IZhang & MacFadvenll200€ 


) to 10" 10 



However, setting this parameter in any problem independent 
way introduces an artificial scale into the problem since dis- 
continuities smaller than this scale are considered to be part 
of the smooth flow. It is therefore preferable to dynamically 
choose e such that it is large enough to avoid the triggering 
of the scheme on machine level noise and yet small enough 
not to in fluence the weights calculation in other cases (see, 
however, IJiang fc Shr3 [l996). We describe the procedure we 
use in Appendix [A3] 

Note that for centre to average and average to centre 
reconstructions some of the optimal weights become neg- 
ativ e, and using n egative weights can lead to an instabil- 
ity (|Shi et al.l l200"2h. This problem can be avoided by a sim- 
ple technique that k eeps the sum of the absolute values of 
the weights bounded l|Shi et al.ll2002T ). We use this technique 
for both average to centre and centre to average reconstruc- 
tions. 

Equation (|24p assigns larger weights to the stencils with 
smaller smoothness indicators, i.e., with flatter reconstruc- 
tion profiles and smaller degree of oscillation. These weights 
become extremely small for a stencil that contains a dis- 
continuity; this is referred to as adaptive stencil choice and 
provides the nonoscillatory property of the scheme. At the 
same time, weights (|24p are designed to be close to the op- 
timal weights according to (|2ip for a smooth flow locally 
well-approximated by a parabola. Therefore, for such flows, 
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the WENO-5 scheme is fifth order in space. However, there 
are other flows for which the scheme becomes third order. 



4.8 Convergence properties of WENO-type 
schemes in smooth flows 

Preserving a high order interpolation near extrema and at 
the same time maintaining the nonos cillatory property o f 
the solution is a challenging problem (|Harten et al.l fl987h . 
For instance, all schemes that satisfy the monotonicity con- 
straint reduce to first order near extrema. 

WENO-5 reconstruction is very appealing in that it 
maintains high order near extrema. For a common type 
of extremum, one with a nonvan ishing second deriv a tive, 
it is c laimed to be fifth order bv Ijiang fc Shul (|l996l ): IShul 
l| 19971 ). However, we have found that in general the WENO-5 
scheme becomes fourth order near such extrema. Further, for 
a smooth flow whose Taylor series expansion is dominated 
by third or higher order terms, the WENO-5 scheme reduces 
to third order. Appendix lA4l discusses the convergence prop- 
erties of the WENO-5 scheme in more detail and gives ex- 
pressions for the smoothness indicators. To help minimize 
the excessive reduction of order by the standard WENO-5 
scheme, our scheme uses the full stencil (|20[) in regions where 
this stencil gives monotonic values for the function and all 
of its derivatives (see Appendix I A5|) . 

One can construct WENO-type schemes of an arbi- 
trarily high order (e.g. seventh, ninth, etc.) that deliver 
an even higher order when the lower-order derivatives van- 
ish. Such schemes provide a much better resolution of 
contact discontinuities and smooth features of the flow 
than thei r lower-order coun terparts, given the same grid 
cell size (jLatini et al.l l2007h . Howe ver, they also re quire 
more computational effort (see, e.g., lOiu fc Shull2002h . Ap- 
pendix IA6I proves some properties of higher order WENO- 
type schemes, discusses their benefits, and points out their 
limitations for handling critical points in smooth flows. 

The use of WENO weights alone is known to be insuf- 
ficient in avoiding spurious oscillations near discontinuities, 
and a known remedy is to use eigenvector decomposition. In 
the next section, we discuss our alternative approach that 
does not require computing the eigenvectors. 



4.9 Oscillation-free reconstruction 

Interpolating the primitive, conserved, or any other arbi- 
trary quantities can in general lead to spurious oscillations. 
Such oscillations can be introduced because the interpola- 
tions allow arbitrary mixing between different eigenmodes. 

One can avoid significant spurious oscillations by recon- 
structing the so-called wavestr engths, the wave amplitudes 
in the local characteristic fields (IQiu fc Shul2002l ). We gener- 
ically refer to this procedure as the eigenvector decomposi- 
tion approach. Interpolating each wavestrength individually 
in most case s leads to only sma ll mixing between different 
eigenmodes l)Harten et al.lfl987"h . There are several types of 
reconstructions based on this idea: characterist ic-wise recon- 
struc t ion, field-by-field reconstruction, e tc. (|Harten et al.l 
Il987l : IShu fc Osherl 1 19881 . 1 19891 : Ishulll997l ). All are compu- 
tationally intensive and require sophisticated coding that 
becomes complicated for the GRMHD case. 



We adopt an alternative approach that involves lo- 
cally red ucing the order of the s c heme near discon- 
tinuities llCole la fc Woodward! Il984l : iLiu fc Osherl Il998l : 
iMignone et al.1 120051 )." The Riemann solver is fed interface 
values that are close to first order near shocks. This proce- 
dure reduces the mixing between different eigenmodes and 
makes the scheme more robust. Further, since this approach 
does not require the computation of any eigenvectors, it can 
be used to solve weakly hyperbolic systems. It is computa- 
tionally more efficient than eigenvector d ecomposition, and 
is easier to implement |Liu fc Osherl 19981 ). See Appendix lA7 
for a discussion. 



5 NUMERICAL TESTS 

We have confirmed the accuracy and robustness of our code 
by running an extensive series of tests. We begin the dis- 
cussion with four tests that show the importance of distin- 
guishing between the average and point values of conserved 
quantities. We then discuss a number of standard tests in 
the literature. Our scheme succe ssfully evolves all nonrela- 
tivistic hydrodynamic tests fro mlLiska fc Wendrofj ([2 003a . 
LW03 ) and all chosen tests from lZhang fc MacFadvenl (|2006 , 
ZM06), but we only discuss those tests we found to be 
most challeng ing. Finally, we discu ss some general relativis- 
tic tests from lGammie et alj l|2003l ). 

Table [1] provides detailed information on some of 
the one-dimensional test problems. Unless ot herwise indi - 
cate d, we use the exact Ri e mann solvers by iTorol (|l997h 
and iGiacomazzo fc Rezzollal (|2006l ) to obtain the exact so- 
lutions for nonrelativistic and relativistic one-dimensional 
Riemann test problems, respectively. 

For nonrelativistic problems we set the value of the 
speed of light in the scheme to 10 10 so that any velocity on 
the order of unity is then nonrelativistic to machine accu- 
racy (i.e. 7 — 1 = 0). All problems use the ideal gas equation 
of state, p g — (r — l)u g , where F is the adiabatic gas index. 

For all tests, we use a constant set of numerical pa- 
rameters that control the behavior of the scheme. This is as 
opposed to many other works that fine-tune their numerical 
parameters in order to make some tests work. The fact that 
we are able to run all the tests with a single set of parameters 
proves the robustness of our method as a single scheme to 
study a vast array of problems. We use WHAM scheme with 
a Courant factor of 0.5 and the approximate HLL Riemann 
solver for all tests. 



5.1 Smooth high Mach number flow: Hubble-type 
flow 

This is a simple high Mach number flow problem that tests 
the ability of a numerical scheme to handle flows with dis- 
parate energy scales. The problem is defined in section [3T2T 
For the test described here the flow parameters are set as 
follows: 



Po = 1, 
v'o = 1, 
wo = 4 X 10~ 



(25) 
(26) 
(27) 
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Table 1. Parameter values selected for the one-dimensional nonrelativistic Riemann problems, shock - entropy wave interaction problem, 
and relativistic Riemann problems. All problems are run on the interval (0, 1), except the shock - entropy wave interaction problem on 
(—5,5), Sod's problem on (—1.5,1.5), and moving Sod's problem on (—1.5,25.5). The initial position of the discontinuity is located at 
x = xq, the tests are run until the final time tp, with gas adiabatic index of T, with N grid cells. Here p is the density, v x and v y the 
3-velocity components, and p g the pressure. 'L' denotes the left state and 'R' the right state. 



Test 


PL 




v y,L 


Pg,L 


PR 




v y,R 


Pg,R 


r 






N 


Sec. 


Sod 


1 








1 


0.2 








0.01 


5/3 


0.0 


0.77 


150 


IB~2l 


Moving Sod 


1 


28.87 





1 


0.2 


28.87 





0.01 


5/3 


0.0 


0.77 


1350 


E3 


Noh 


1 


1 








1 


-1 








5/3 


0.5 


1 


100 


rp eh 
15.51 


Moving contact 


1.4 


0.1 





1 


1 





0.1 


1 


1.4 


0.5 


2 


100 


IS~6l 


Test 4 


5.99924 


19.5975 





460.894 


5.99242 


-6.19633 





46.095 


1.4 


0.4 


0.035 


200 


IS~7l 


Shock - entropy 


3.857143 


2.629369 





10.33333 


l + 0.2sin(5x) 








1 


1.4 


-4 


1.8 


400 


15^1 


Rel. problem 1 


10.0 








13.33 


1.0 








lO" 8 


5/3 


0.5 


0.4 


400 


|5.14| 


Rel. problem 2 


1.0 








1000 


1.0 








lO" 2 


5/3 


0.5 


0.4 


400 


|5.15| 


Rel. problem 3 


1.0 


0.9 





1 


1.0 








10 


4/3 


0.5 


0.4 


400 


|5.16| 


Rel. problem 4 


1.0 








1000 


1.0 





0.99 


10" 2 


5/3 


0.5 


0.4 


400 


|5.17| 


Rel. problem 5 


1.0 





0.9 


1000 


1.0 





0.9 


10" 2 


5/3 


0.5 


0.6 


400 


|5.18| 


Rel. problem 6 


1.0 


1-10" 10 





0.001 


(uses reflecting boundary at 


x = 1) 


4/3 


1.0 


2 


100 


|5.19| 



with r = 1.4. We employ a resolution of 64 cells on the in- 
terval (—1, 1), so that the Mach number varies from M m i n ~ 
104.4 to M max « 6577.1 at the grid cell centres. For the val- 
ues in the boundary cells we use linearly extrapolated values 
of the primitive quantities. 

The Hubble-type flow problem illustrates how crucial 
it is to account for the difference between cell averaged and 
cell centred conserved quantities in high Mach number flows. 
Table [2] shows the relative Li-error in the internal energy 
and other quantities at the characteristic time of evolution, 
t F = 1, for the WHAM, WENO-IFV, HARM0 and Athenffj 
schemes. We define the absolute Li-error norm of any quan- 
tity it as 

Aau ee (Au) Ll = ]T | u — ical - u f* ct l/AT, (28) 
i 

where N is the number of elements, and the relative Li-error 
norm as 

A R u = (Att/it max )L 1 = A A ?V max l u r act |- (29) 

i 

Since WHAM performs the proper conversion between cell 
averaged and cell centred conserved quantities before the 
reconstruction step, it is 4th order accurate in time and gets 
the final distribution of internal energy with a relative Li- 
enor of only 5% at the default resolution of A = 640 

In contrast, the WENO-IFV, HARM, and Athena 

4 The HARM scheme converges at second order in space and 
time and does not account for the difference betwe en cell averaged 
and ce ll centred values of conserved quantities, see lGammie et aU 
ll2003r i. 

5 We use Athena 2.0 ( Stonel l2006l) with third order spatial inter- 
polation and the Roe solver. We use a Courant factor of 0.5. For 
setting the values in the boundary grid cells we use quadratically 
extrapolated values of conserved quantities. 

6 4th order accuracy in time at one time step implies an error 
term of 0(At 5 ). Therefore, the error due to evolution from t = 
to t = t F is n X 0(At 5 ) ~ 0(Ai 4 ) ~ 0(C 4 ), where n = t F /At is 
the number of time steps to reach the final time tp and C = O(At) 
is the Courant factor. 



Table 2. Relative L\ -error and convergence order in inter- 
nal energy at the final time tp = 1 for the Hubble-type flow 
problem (section 15.111 . The Courant f actor is C, and the reso- 
lution is N. The WE NO-IFV, HARM llGammie et al.ll2003h . and 
Athena llStonell2006T) schemes converge at second order in space 
and at zeroth order in time in agreement with equation 113t . 
One has to increase the resolution by a factor of more than 100 
to match the performance of the WHAM scheme. WHAM con- 
verges at 4th order in space and time, see the note after equation 

eg. 



Scheme C N Aru 9 Order 



WHAM 


0.5 


m 


|5.2e-2| 






0.5 


128 


4.3e-3 


3.6 




0.5 


256 


3.1e-4 


3.8 




0.5 


512 


2.1e-5 


3.9 




0.5 


1024 


1.4e-6 


3.9 




0.05 


64 


5.2e-6 


4.0 


WENO-IFV 


0.5 


64 


1278 






0.5 


640 


12.78 


2.0 




0.5 


6400 


12.78c-2 


2.0 




0.5 


| 9600 | 


|5.7c-2| 


2.0 




0.5 


25600 


0.8c-2 


2.0 




0.05 


64 


1278 


0.0 


HARM 


0.5 


64 


773.3 






0.5 


640 


7.8 


2.0 




0.5 


6400 


7.8c-2 


2.0 




0.5 


| 7680 | 


|5.4e-2| 


2.0 




0.5 


9600 


3.4e-2 


2.0 




0.05 


64 


1273 


-0.2 


Athena 


0.5 


64 


2527 






0.5 


640 


30.43 


1.9 




0.5 


6400 


0.14 


2.3 




0.5 


| 12000 | 


|4.2c-2| 


1.9 




0.05 


64 


2459 


0.01 
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-5xl0 2 



0.2 0.4 0.6 

(b) HARM 



0.8 



Figure 3. Analytical and numerical solutions for the Hubblc- 
type flow problem, section 15.11 at a resolution of N = 64. At 
each time step a scatter plot of dimensionless values of internal 
energy is shown for WHAM (a) and HARM (b). The analytical 
solution H 12 ^ for internal energy is shown with a solid line on each 
panel. Note that on panel (b) the analytic solution is compressed 
into a horizontal line because the vertical scale has been changed 
by a factor of about 1000. WHAM accurately captures the de- 
pendence of the internal energy in this highly supersonic flow, 
while other schemes (including HARM and Athena) that do not 
differentiate between the average and point values of conserved 
quantities makes a large error in internal energy. 



codes, which do not perform the average to centre conver- 
sion, produce large errors. In fact, they obtain a negative 
value of the internal energy already after the first time stepQ 
When the internal energy is negative, for the WENO-IFV 
and HARM schemes we set the sound speed to 0; for Athena, 
we use its default settings. 

Figure [3] visualizes the time dependence of the internal 
energy for the WHAM and HARM schemes. The solution for 
internal energy given by WHAM deviates only slightly from 
the analytic solution. The internal energy given by HARM 
becomes negative after the first time step and is very nonuni- 
form in space. Any scheme (including HARM and Athena) 
that does not differentiate between the average and point 
values of conserved quantities will make a similar error. This 
error in the internal energy is second order in space, with 
a large coefficient proportional to the square of the Mach 
number (see eq. I13|l . This error does not decrease if the time 
step is decreased. For these schemes one has to use a reso- 
lution more than 100 times larger than the default one to 
match the accuracy of the WHAM scheme at the default 



7 We note that a first order Godunov scheme with an approxi- 
mate HLL Riemann solver produces inaccurate but physical val- 
ues of the internal energy. 




(b) 



Figure 4. Panel (a) shows the standard Sod's problem. Panel (b) 
shows the 'moving Sod's problem' with the pre-shock state mov- 
ing supersonically through the grid at a Mach number of 100. 
The numerical solution by the WHAM scheme is shown with con- 
nected dots and the analytic solution is shown with a solid line. 
For the lower panel (b) the ^-coordinate has been remapped to 
correspond to the same range as in panel (a). WHAM is able to 
accurately treat shocks when both the pre-shock and post-shock 
regions are moving supersonically. 



resolution (see table [2]). In other words, for high Mach num- 
ber problems, high order schemes are much more effective 
than lower order ones. 



5.2 Discontinuous high Mach number flow: Sod's 
shock tube 

In t his section we co nsider two versions of Sod's prob- 
lem (|Trac fc Penll2004f ). The first 'stationary' version is a 
standard, simple shock tube problem, a variation of the orig- 
inal problem by ISodl (1 19781 ) . Table [1] shows the initial con- 



ditions for this problem and figure 4(a) shows the density 



distribution given by the WHAM scheme overplotted on the 
exact solution at the final time. 

The second 'moving' version is the 'stationary' problem 
boosted to a supersonic speed such that the pre-shock Mach 
number of the initial right state of the problem is equal to 
100. This verifies the code's ability to handle shocks with 
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both pre-shocked and post-shocked gas moving supersoni- 
cally with respect to the grid. Godunov-type schemes that 
operate on a fixed Eulerian grid are claimed to have serious 
difficulties, e.g., producing the wrong height of the shock 
and g enerating spurious post-shock oscillations (|Trac fc Pen! 
12004 ). 

In studying the 'moving' version of Sod's problem with 
WHAM, we keep the same grid cell spacing and the final 
time as in the 'stationary' version but extend the grid so 
that the wave structure does not go off the grid until the final 
time of the test (see table [TJ. We have performed this prob- 
lem wit h both the WHAM sc heme, see figure 4(b) and with 
HARM (|Gammie et al.ll2003h . We find that neither scheme 
sh ows the sugg ested violent post-shock oscillations claimed 
bv lTrac fc Peru . It could be that their numerical scheme has 
difficulties because it reconstructs the conserved quantities, 
not primitive quantities. The lower order HARM scheme 
gets an incorrect height of the shock for the 'moving' prob- 
lem (and produces the correct shock height in the 'station- 
ary' version of the problem) . This phenomenon also occurs 
with the Eulerian scheme of iTrac fe Per] (12004 ) 

W ithout the moving grid technique of ITrac fc Per] 
|2004h . WHAM is able to obtain a comparable resolution 
of the contact discontinuity and value of the shock height. 
This shows that accounting for the difference between the 
point and average values of conserved quantities is an a l- 
ternative to the moving grid scheme of lTrac fc Per] (2004). 
WHAM converges uniformly at first order for this test. 

5.3 One-dimensional hydrodynamic caustics 

In this problem (see figure [5} a highly supersonic smooth 
flow steepens into a pair of shocks posing challenges for nu- 
merical schemes: the ability to accurately evolve a smooth 
high Mach number flow and handle its interaction with 
strong shocks. Initially the density and the pressure are uni- 
form, p = 1 and p = 10 -10 , and the velocity is a sine wave, 
v(x) = (27r) _1 sin(27rj:). This corresponds to Mach number 
reaching 1.2 x 10 4 and the internal energy accounting only 
for about 10~ 8 of the kinetic energy. We run the problem on 
an interval (0, 1) until time tp = 3 with periodic boundary 
condi tions and we u se a resolution of 128 grid cells as chosen 
bv lRvu et al] (| 19931 ). 

For this problem iRvu et all (|l993l ) find that schemes 
evolving the total energy make unacceptably lar ge errors in 
the in ternal energy even at very high resolutions. I Ryu et all 
(|l993h therefore argue for the 'dual-energy formalism' in 
which one switches to evolving the entropy equation in the 
high Mach number regions of the flow. Although such an ap- 
proach does not suffer from the high Mach number problem 
in smooth supersonic flows, it cannot account for dissipation 
in weak shocks (see section [TJ. 

WHAM correctly captures shocks in highly supersonic 
flows (see the previous section) and for smooth flows quickly 
converges to the correct solution as the resolution is in- 
creased. Figure [5] shows WHAM getting the final distribu- 
tion of pressure in the high Mach number region within a 
factor of 3 from the converged solution. This is 4 orders 
of magnitude be tter than the energy-conserving scheme of 
IRvu et all (|l993l ). The shock is always resolved with about 
four points. Using a resolution of 180 and a Courant factor 
of 0.2 gives a pre-shock pressure error of only 10%. With 



a resolution of 256 and the standard Courant factor of 0.5, 
the error is reduced by a factor of 10 ~ 2 3,5 to 25%, i.e. our 
scheme converges in pointwise sense at 3.5th order in the 
pre-shock region. HARM or Athena, which are strictly sec- 
ond order in supersonic flows and do not improve their error 
for smaller Courant factors, obtain an error ~ 25% only at 
a resolution of ~ 45, 000 grid cells. 



5.4 Linear wave advection 

To verify the order of convergence of WHAM, we perform 
the advection of smooth sound and density waves in Carte- 
sian coordinates in two dimensions. This verifies the order 
and accuracy of the multidimensional reconstruction and of 
the Runge-Kutta time stepping. 

We have set up this test in the same way as in the 
Athena code |Stonell2006l ). For both the density and sound 
wave advection problems, we choose a wave with a wave 
vector k = (k x ,k y ) = (2tt/ sin a, 2n/ cos a) directed at an 
angle a — tan _1 (2) « 64.4 degrees w.r.t. the a;-axis. We use 
a computational box of size (0,sina) x (0,cosq) such that 
periodic boundary conditions can be applied. The number of 
grid cells in the i-direction is twice that in the y-direction so 
that each grid cell is a square even though the computational 
box is a rectangle. Therefore the wave does not travel along 
the diagonal of grid cells guaranteeing that the test is truly 
multidimensional. 

We use an equation of state with F = 5/3 and choose a 
uniform background state with a unit density po = 1 and a 
unit sound speed c s ,o = 1, which corresponds to the back- 
ground internal energy no = 0.9. For the sound wave test we 
choose the background velocity to be zero (vo = 0) and for 
the density wave test we choose the background velocity to 
be unity along the direction of the wave vector (vo = k/|k|). 
We run both tests for one temporal period T of the wave, 
i.e. until tp =T = 0.4. Relative perturbations for the sound 
wave include perturbations in density, velocity, and internal 
energy, 



Sp/po 
5v/c 8 , 
Su/uo 



= Acos (k • r) 




(30) 



whereas for the density wave the perturbations are only in 
density, 



Sp/po 
5v/c 8 , 
Su/uo 



= Acos (k • r) 




(31) 



We choose the relative magnitude A of the perturbation such 
that nonlinear wave effects for the sound wave at the final 
time lead to relative nonlinear corrections 5nl to the solution 
on the order of the numerical precision for our machine, 
<5nl ~ A 2 t F /T = e machinc » 10" 15 , so we set A = 3.2 x 10" 8 . 

We perform the calculations at different resolutions. Ta- 
ble [3] shows that the WHAM scheme indeed converges at 
fifth order. For comparison, we have also run the density and 
sound wave tests using the WENO-IFV scheme, HARM, and 
Athena. These three schemes converge only at second order 
despite the cell centre to interface reconstruction used in the 
WENO-IFV scheme being fifth order. This shows that not 
making a distinction between the average and point values 
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Figure 5. Density, pressure, and velocity for the one-dimensional hydrodynamic caustics problem at time tp = 3. Connected dots show 
the result of WHAM at the resolution of N = 128 grid cells with the standard Courant factor of C = 0.5 and the solid line shows the 
converged solution obtained with N = 2048 grid cells. With N = 128 and C = 0.5, WHAM reproduces the low pressure pre-shock region 
with an error of 250% that is several orders of magnitude smaller than with the energy conserving scheme of iRvu et all il993h . With 
N = 180 and C = 0.2, the pre-shock pressure error is only 10% with WHAM, and with N = 256 and C = 0.5 the error is 25%. Highly 
accurate, higher-order energy conserving schemes like WHAM are capable of obtaining accepta ble solutions for extremely high Mach 
number flows without resorting to the dual-energy formalism, recommended bv lRvu et al.l [|l993T ) 



Table 3. Relative Li-error in density, A^jp, and the order of 
convergence (base 2 logarithm of the error) of WHAM, WENO- 
IFV and HARM schemes for the two-dimensional density and 
sound wave problems (section l5,4[ l. The WHAM scheme converges 
to the analytic solution at fifth order for both problems, while 
the WENO-IFV scheme and the HARM scheme converge only at 
second order. 



Resolution 


16 x 8 


32 x 16 


64 x 32 


128 x 64 




Entropy wave 






WHAM 


3.9e-09 


9.2&-11 


2.6e-12 


8.0e-14 


Order 




5.5 


5.2 


5.0 


WENO-IFV 


5.5e-09 


7.1c-10 


1.7e-10 


4.4c-ll 


Order 




3.0 


2.0 


2.0 


HARM 


7.7e-09 


1.6e-09 


5.9e-10 


1.7c-10 


Order 




2.2 


1.5 


1.7 


Athena 


2.9e-09 


l.lc-09 


2.7e-10 


6.9c-ll 


Order 




1.4 


2.1 


2.0 




Sound wave 






WHAM 


4.0e-09 


9.2e-ll 


2.6e-12 


8.0e-14 


Order 




5.4 


5.1 


5.0 


WENO-IFV 


5.5e-09 


7.1e-10 


1.7e-10 


4.4c-ll 


Order 




3.0 


2.0 


2.0 


HARM 


7.9e-09 


1.8c-09 


6.3e-10 


1.9e-10 


Order 




2.2 


1.5 


1.8 


Athena 


2.5e-09 


6.5e-10 


1.8e-10 


4.2c-ll 


Order 




1.9 


1.9 


2.1 



of conserved quantities and fluxes can significantly impact 
the performance of schemes even for problems at a low Mach 
number. We have als o perfo rmed a smooth density advec- 
tion test described by ILW03I and found similar convergence 
results. 




Figure 6. One-dimensional Noh problem. The analytic solution is 
shown with a solid line, and the numerical solution is shown with 
connected dots. WHAM does not exhibit any Gibbs phenomenon 
near the two strong shocks and has a very small dip in density at 
the centre. 



5.5 Noh 

This is a one-dimensional Riemann problem that tests the 
ability of codes to handle infinite-strength shocks and tests 
the faithfulness of reproduction of shock-shock interaction. 
Initially, two streams are plunging into each other symmet- 
rically with a constant velocity. As a result of their interac- 
tion, two shocks develop that travel away from each other 
and leave matter behind at rest with a constant density and 
pressure (see figure [5} . This test problem shows the impor- 
tance of reducing to lower order in strong shocks. If we run 
this problem without any stencil reduction or eigenvector de- 
composition (see section |X3J) , the result exhibits significant 
Gibbs phenomenon. Even the second order HARM scheme, 
which uses the MC limiter for spatial reconstruction, pro- 
duces post-shock oscillations in density with an amplitude 
of about 5%. With WHAM, such post-shock oscillations are 
absent and the dip in density at the centre is also small. 
Our scheme's result is comparable to that of the WENO- 
5 scheme and is superior to that of PPM from the study 
of lLWOi The WENO-5 scheme uses field-by-field decompo- 
sition that helps it avoid oscillations (see section T3.4p ; the 
PPM sche me uses a form of stencil redu ction (so-called flat- 
tening, see IColella fc Woodward! \l984) that locally lowers 
the order of the scheme to first order near shocks. For the 
latter scheme the dip at the centre is large and some post- 
shock oscillations are visible. 
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Figure 7. Moving contact problem. The analytic solution is 
shown with a solid line, and the numerical solution is shown with 
connected dots. WHAM exhibits only a moderate diffusion of the 
contact discontinuity despite avoiding eigenvector decomposition. 



5.6 Moving contact 

Figure [7] shows a test that measures the amount of dif- 
fusion of contact discontinuities in the numerical scheme. 
In the exact solution the initially sharp contact disconti- 
nuity remains sharp throughout the evolution. Used alone, 
the component-wise WENO-type reconstruction would lead 
to excessive smearing of the contact discontinuity because 
it does not use the full stencil inside the steeply changing 
density profile: it chooses the one-sided stencil inside of the 
discontinuity. To avoid this excessive smearing, we force the 
use of the full stencil in the regions where the polynomial 
fit due to the full stencil and all of its derivatives are mono- 
tonic (see Appendix I A5|) . Our result is comparable to that 
of W ENO-5 scheme, which uses eigenvector decomposition, 
from lLWOl 




Figure 8. Test 4 from |lW03. The analytic solution is shown 
with a solid line, and the numerical solution is shown with con- 
nected dots. The figure shows that WHAM accurately resolves the 
structure of this complicated Riemann problem. In particular, it 
provides good resolution of the moving contact discontinuity. 




Figure 9. Snapshot of the density distribution for the shock - 
entropy wave interaction problem at the final time. The converged 
solution (at the resolution of 2000 grid cells) is shown with the 
solid line, and the numerical solution is shown with connected 
dots. WHAM accurately resolves the interaction of the shock and 
the smooth flow. 



5.7 Test 4 problem from ILW03I 

In this problem the high resolution (small width) of the 
contact discontinuity, located in figure [8] at x ~ 0.7, and 
the absence of oscillations in the two high density regions 
are important. The numerical solution is complicated by 
the fact that both of these regions are 'built-up' as a re- 
sult of the decay of an initial discontinuity. The HARM 
scheme exhibits post-shock oscillations at the level of about 
5% in the left built-up state. WHAM does not exhibit such 
Gibbs phenomenon. The resolution of the contact discon- 
tinuity in ou r schem e is as good as that of the WENO-5 
scheme from LW03, which requires eigenvector decomposi- 
tion. The resolution of the contact by WHAM is the same or 
better than that of other schemes studied bv lLW03t except 
the PPM scheme, which u ses an artificial con tact discontinu- 
ity sharpening technique l|Frvxell et al.ll2000|). a Lagran gian- 
remap version of PPM called VH1 (jBlondin et al.lll99ll ). and 
the WAFT scheme that uses the HLLC Riemann solver. 



comp arable result t o that of the CWENO- 5 l|Qiu fc Shul 
12002 1 and WENO-5 (jLiska fc Wendrofi1l2003bl ) schemes that 
use this decomposition. 



5.9 Nonrelativistic 2D Riemann problem 

We have run all two-dimensional test problems from LW03 
and the results of our code are comparable to those of other 
codes presented there. We picked the particular test prob- 
lem 4 to show here because it is the only one that exhibited 
any noticeable Gibbs phenomenon. The initial conditions 
for this test problem are shown in table [4] The problem is 
computed at a resolution of 400 x 400 and uses T = 1.4. 
This test problem initially contains 4 planar shocks. Fig- 
ure [10] shows the final state of the problem where a high- 
density eye-shaped area develops, bounded by two shocks. 
Even though stationary contacts are slightly less resolved in 
our code than in other schemes, for moving types of contact 
we are the same or more accurate. 



5.8 Shock — entropy wave interaction test problem 

This test (figure [9} challenges the ability of the code to han- 
dle the interaction of a shock and a smoothly varying flow. 
Most second order schemes (jLiska fc Wendrofll l2003bl ) in- 
cluding HARM provide inadequate resolution of the wave 
structure resulting from the interaction of the shock with the 
stationary density wave. WHAM is able to accurately resolve 
the high-frequency waves that develop behind the shock and 
provides, without the use of eigenvector decomposition, a 



5.10 2d Noh 

This is a 2-dimensional version of the Noh problem for an 
ideal gas with F = 5/3 set up in Cartesian coordinates. 
Initially, the flow is cylindrically symmetric and converging 
on to the origin with a constant radial velocity, v r = 1. The 
density distribution is uniform, p — 1, and the pressure is 
zero. Since the flow converges to a point, an outgoing shock 
develops as can be seen in figure 1111 The shock position is 
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Table 4. Initial conditions of the n onrelat ivistic two-dimensional 
Riemann problem, test case 4 from lLWOl The final time for this 
test problem is tp = 0.25, and T = 1.4. The upper row of the 
table lists the initial state of the upper left and right corners of 
the Riemann problem: the 'L' and the 'R' indices stand for the 
left and right states, respectively. The lower row lists the lower 
two states. Otherwise the notation is the same as in table [l] 



PL 


V X ,L 


V V,L 


PL 


PR 


v x ,R 


v y,R 


PR 


0.5065 


0.8939 


0.0 


0.35 


1.1 


0.0 


0.0 


1.1 


1.1 


0.8939 


0.8939 


1.1 


0.5065 


0.0 


0.8939 


0.35 




Figu re 10 . 2D nonrelativistic Riemann problem, case 4 
from lLW03l . see section 15.91 Pressure is shown in color (red de- 
notes high values and blue low values) overplotted by a set of 
density con tours g oing from 0.52 to 1.92 with a step of 0.05, the 
same as in LWOjJ for easier comparison. WHAM accurately re- 
solves the density structure in the eye-shaped area. The minor 
oscillations are comparable to other published results. 



located at R = R a = t/3, where R = \/ x 2 + y 2 . In the post- 
shock region (R < R s ) the density is p — 16, the velocity is 
zero and the pressure is p g — 16/3. In the pre-shock region 
(R > Rs) the density is p — 1 + t/R, the vel ocity stays 
constant |v| = 1, and the pressure remains zero l|Noh|| 19871 ; 
iLiska fc Wendrofj|2OT3ah . The initial conditions are evolved 
until tp = 2. We use a resolu tion of 400 x 400 cells. Further, 
for consistency with iLWOSl we initialize the pressure with 
a small value p g — po — 10~ 6 so that it is dynamically 
unimportant. 

This test problem is unusual in that most of the final 
state is determined by the time-dependent boundary condi- 
tions. In contrast to the one-dimensional Noh problem, here 
the boundary conditions have to be the evolved initial con- 
ditions. For each of the Runge-Kutta substeps we set the 
time we use for the boundary conditions to correspond to 
the time of that trial time step. We find that using p — po 
for the boundary condition as is done in, e.g.. |LW03| , results 



in a sharp kink in pressure along the diagonal in the region of 
influence of the boundary conditions. To avoid this, instead 
of keeping the pressure constant at the boundary, we use an 
approximate solution for the pressure p g = po(l + Tt/R) in 
the pre-shock region. We obtain this approximate solution 
by solving the internal energy advection equation, 

with an initial condition p g = po and with all other quan- 
tities determined by the above analytic solution (i.e. we ne- 
glect the effect of the pressure force on the evolution of den- 
sity and velocity). We have not found a closed analytic so- 
lution to this problem for a nonzero initial pressure. 

The smooth pre-shocked region is a highly supersonic 
flow with an initial Mach number of M « 775. For the evolu- 
tion of the internal energy to be accurate in this supersonic 
region, the truncation error has to be small. We find that for 
this both the de-averaging of the conserved quantities and 
the transversal averaging of the fluxes are important. If we 
do not perform either of these two operations, the internal 
energy in the smooth region becomes negativ e. We expect 
that any scheme, such as PPM used in ILWOol . that ignores 
the difference between points and averages will generate a 
significant error in the internal energy in the pre-shock re- 
gion. 

The numerical results are shown in figure [TT] and are 
much superio r to th ose of HARM and the other schemes 
considered bv lLW03t The numerical solution is very smooth 
in the pre-shock region and does not show any visible os- 
cillations for radii far enough (r > 0.2) from the origin in 
the post-shock region except near the shock. The solution 
for density remains accurate even if we use p g = po for the 
boundary condition. This should be contrasted to the results 
of other codes, all of which show much more significant os- 
cillations both in the pre-shock and post-shock regions. 

5.11 Implosion 

The implosion test problem (iLWOot ) corresponds to the 
interaction of a low-density, low-pressure diamond {pi = 
0.125, pi = 0.14) with a high-density, high-pressure exterior 
{po — 1, Po = 1) in a square box with reflecting bound- 
aries, both the diamond and the box centred on the ori- 
gin. The gas adiabatic index is V = 1.4 The vertices of 
the diamond are located at the intersections of coordinate 
axes with a circle of radius 0.15. The computational box is 
(—0.3, 0.3) x (—0.3, 0.3). Initially the velocities are zero. We 
perform the co mputati ons only for the upper-right quadrant 
of the box (as in lLW03l we use reflecting boundary conditions 
on all 4 boundaries) at a resolution of 400 x 400. In order 
to avoid grid-induced artifacts in the initial conditions, for 
grid cells intersected by the discontinuity we use the average 
of the two states. The initial discontinuity evolves into two 
shocks with a contact discontinuity between them. 

The interaction of the contact discontinuity with the 
shocks generates streams that travel toward the origin and, 
after colliding there, form an outflow in the form of a nar- 
row jet. Figure [12] shows the configuration at the final time 
tp = 2.5. This problem tests the ability of the code to re- 
solve contact discontinuities and maintain symmetry in two 
dimensions: if a code does not preserve symmetry, the jet will 
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Figure 11. Density distribution in the two-dimensional Noh problem at the final time tp = 2. The left panel shows density both in color 
(red denotes high val ues and blue low values) and by a set of contours going from 2.5 to 4 in steps of 0.25 and from 14 to 17 in steps 
of 0.2, the same as in lLWOij . The right panel shows the scatter plot of density w.r.t. radius (blue dots) and the analytic solution (thin 
red line). WHAM provides a very accurate solution in the pre-shock region and shows significantly fewer oscillations behind the shock 
compared to other schemes. For the left panel, the noise near the shock appears because one of the contours is chosen to the analytic 
value. 



not be produced or will be distorted. Lower order schemes 
significantly diffus e the je t at this resolution. Unlike some of 
the schemes from HA/VPS, our code gives exactly symmetric 
results. Further, the resolution of the jet provided by our 
code i s comparable to that of the WENO-5 scheme from 
and is superior to all other schemes discussed in that 
paper and HARM. Since Athena has a lower dissipatio n than 
WHA M, their jet head travels about 10% further (IStond 
l200fj ). The result obtained with WHAM could be improved 
if we use an approximate HLLC Riemann solver, which pro- 
vides a better resolution for contact discontinuities. 



5.12 Explosion 

The explosion problem verifies the ability of the code to 
evolve unstable contacts. In particular this problem studies 
how sensitive the code is to numerical perturbations, which 
arise from the discreteness of the grid. 

The initial conditions are cylindrically symmetric, with 
a high-density, high-pressure cylinder, pi = 1, pi = 1, with 
a radius of 0.4, embedded in a low-density, low-pressure 
medi um, p = 0.125, p = 0.1. The gas constant is T — 1.4. 
As in ILWOSl we perform the computations in Cartesian co- 
ordinates in a square box (0, 1.5) x (0, 1.5) at a resolution of 
400 x 400 grid cells, with reflecting boundary conditions at 
the left and bottom boundaries and zero-derivative condi- 
tions at the upper and right boundaries. The initial discon- 
tinuity evolves into two shocks with a contact discontinuity 
in between. Figure [13] shows a snapshot of the problem at 
a final time of tp = 3.2 after the outgoing shock leaves the 
computational domain through the upper and right bound- 
aries and the ingoing shock bounces off the origin and passes 
through the contact discontinuity. 




Figure 12. Implosion problem l|LW03h . Pressure is shown in color 
(red denotes high values and blue low values) overplotted by a 
set of density contour s going from 0.35 to 1.1 in steps of 0.025, 
the same as in LW03. The figure shows that WHAM is able to 
maintain perfect symmetry of the solution and has low diffusivity 
that is required for producing the narrow jet. 



The interface of the contact discontinuity in this test is 
unstable to perturbations, so the initial conditions have to 
be properly averaged for grid cells that are intersected by 
the discontinuity to avoid seed perturbations in the contact 
discontinuity coming from the grid structure in the initial 
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Figure 13. Explosion problem (see, e.g., llW03h . Pressure is 
shown in color (red denotes high values and blue low values) 
overplotted by a set of density c ontou rs going from 0.08 to 0.21 
in steps of 0.005, the same as in ILW03I . WHAM shows the same 
or smaller degree of breakup of the unst able contact discontinuity 
compared to the WENO-5 scheme from lLW03l . 

conditions. However, one cannot fully avoid perturbations 
to the contact discontinuity, and the discontinuity breaks 
up simil arly to the way it ha ppens fo r other high order 
schemes (jlW03l ). As pointed by[LW03, this problem is sen- 
sitive to the conditions that are applied at the upper and 
right boundaries. Therefore, the interaction of these bound- 
ary conditions with the flow and, in particular, the contact 
discontinuity could seed additional perturbations to the con- 
tact discontinuity. The width and degree of break up of the 
contact discontinuity is the sa me as or smaller than that of 
the WENO-5 scheme in lLW03l . However, unlike the WENO- 
5 scheme, our result is obtained without the use of eigenvec- 
tor decomposition. 



5.13 Moving Gresho problem 

This problem tests how well a numerical scheme is able to 
advect a smooth vortex supported by the balance of pres- 
sure and rotation. A non-moving vortex centred at (0, 0) has 
the following distribution of velocity and pressure in polar 
coordinates (r, <p) with the origin at its centre: 

v v = 5r, p g =5 + 25 / 2 r 2 , 0.0 sC r < 0.2, 

v v = 2-5r, p g = 9 + 25 / 2 r 2 - 20r + 4 In 5r, 0.2 sC r < 0.4, 
%=0, p 9 = 3 + 41n2, 0.4 sCr, 

the radial velocity v r is zero, the density is unity, and the 
polytropic index of the gas is Y — 1.4. In the test, this vor- 
tex is initially imparted unit velocity in the rE-direction. Fig- 
ure [14] shows snapshots of the vortex at the initial time and 
the final time tp = 3 when the vortex centre has moved to 
(3, 0). We perform the computation in cartesian coordinates 
in the rectangle (—0.5,3.5) x (—0.5,0.5) at a resolution of 



Table 5. Relative L\ -errors for the Moving Gresho problem 
(sec. I5.13t . The errors are computed over the core of the vortex 
within radius 0.2 from the vortex centre, at the final time tp = 3. 
WHAM makes about a factor of 5 smaller error in baryon density 
and internal energy than HARM at the resolution of 160 X 40. 



Scheme 






ArV x 


A R Vy 


WHAM 


1.68e-03 


3.90e-03 


2.10e-02 


1.09e-02 


HARM 


8.36e-03 


2.61e-02 


3.87e-02 


2.06e-02 




(b) Final WHAM (c) Final HARM 



Figure 14. Moving Gresho problem (see, e,g.. llW03h . Initial (a) 
and final (b) & (c) distributions of vorticity are shown in color 
(red denotes high values and blue low values) overplotted by a 
set of density c ontour s going from 0.97 to 1.03 in steps of 0.006, 
the same as in LW03. The resolution of the fragments shown is 
40 X 40, the resolution of the full grid is 160 X 40. The WHAM 
scheme (panel b) preserves the structure of the vortex in density 
and velocity exceptionally well co mpared to HARM (panel c) and 
the various schemes considered in llWOj . 



160 x 40 and use zero-derivative outflow bound ary con ditions 
at all boundaries. This is the same as used by[LW03. 

WHAM does extremely well in preserving the structure 
of the vortex in both velocity and density. Our code main- 
tains the vorticity and makes a very small error in density 
so that it ends up with much fewer densit y contours com- 
pared to HARM and all sc hemes in llW03l . HARM as well 
as the codes considered in llW03l such as WENO-5, PPM, 
VH1, etc., destroy the shape of the vortex, making it ellip- 
soidal and/or significantly diffuse in the radial profile. This 
shows the superiority of our high accuracy finite-volume ap- 
proach in applications involving smooth problems. Relative 
Li-errors of the final solution at £f = 3 w.r.t. to the analytic 
solution for WHAM and HARM are shown in table [5] 
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Figure 15. The distribution of density at the final time o f the 
one-dimensional relativistic Riemann problem 1 from ZMQfJ The 
numerical solution is shown with connected dots, and the ana- 
lytical solution is shown with the solid line. WHAM correctly re- 
produces all ingredients of flow: the rarefaction wave, the contact 
discontinuity, and the shock. 




Figure 16. Density distribution for the one-dimensional rela- 
tivistic Riemann problem 2 from ZM06. The numerical solu- 
tion is shown with connected dots, and the analytical solution 
is shown with the solid line. WHAM performs well by getting the 
height of the built-up state to about 70% of the analytic value , 
a height compa rable to other codes l|Zhang fc MacFadver]|2006t 
IMorsonv et~al]|2006h . 



5.14 ID relativistic Riemann problem 1 

The initial conditions are given in table [1] This problem ver- 
ifies the ability of a numerical scheme to treat basic relativis- 
tic problems. As seen in figure [151 the initial discontinuity 
decays into a rarefaction wave, a contact discontinuity, and 
a shock. The resolution of the contact discontinuity and the 
amount of post-shock oscillat ion is comparable to that of 
the WENO scheme in ZM06. Note that our code achieves 
a comparable result despite not using the characteristic de- 
composition that certainly helps the resolution of Riemann 
problems. Table [6] shows the Li-errors of the solution and 
the order of convergence of WHAM for this problem as well 
as the following five relativistic Riemann problems. 



Table 6. Absolute and relative L\ errors in density and the order 
of convergence of WHAM for relativistic Riemann problems 1-6. 
For a unit interval, at which these problems are set up, the defi- 
nition of the absolute L\ error (1 28 I t is equivalent to the definition 
used by IZMOd WHAM converges at approximately first order 
for all tests, where the errors and convergence rates are similar 
to other published work. 
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5.15 ID relativistic Riemann problem 2 

The data for this test are given in table [U This is a chal- 
lenging test with sparse resolution. It is very hard to get the 
correct height of the built-up state due to the smearing of 
the contact discontinuity and the proximity of the shock (fig- 
ure!16ll. Comparable to that of other codes ( Alov et alii 19991 ; 
IZhang fc MacFadvenl 120061 : IMorsonv et all l2006t l. WHAM 
gets the height of the built-up state to about 70% of the 
analytic value. 



5.16 ID relativistic Riemann problem 3 

The initial discontinuity in this problem generates two 
shocks and a contact discontinuity, see figure 1171 Table Q] 
lists the problem parameters. The slow moving reverse shock 
appears to be hard for most codes to handle. Even the most 
dissipative second order schemes like MINMOD produce os- 
cillations at a noticeable level. Our code exhibits them too 
at about twice the level of second order codes that do not 
use eigenvector decomposition. We note that the relativistic 
version of the FLASH code, which uses a flattening pro- 



© 2006 RAS, MNRAS OOO.ITU331 



20 A. Tchekhovskoy, J. C. McKinney, & R. Narayan 



6 


II 1 1 II 1 1 II 
1 1 1 




1 


III 

1 


1 


P 4 












2 














: 


i i 


1 




r-H : 



0.2 0.4 0.6 0.8 1 
x 



Figure 17. Density distributio n of th e one-dimensional relativis- 
tic Riemann problem 3 from IZMOESl . The numerical solution is 
shown with connected dots, and the analytical solution is shown 
with the solid line. This figure shows that WHAM gets the po- 
sitions of discontinuities correctly but sometimes produces small 
amplitude oscillations near slowly moving shocks. Similar oscil- 
lations are present at some level in most codes that do not use 
eigenvector decomposition 
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Figure 18. One-dimensional relativistic Riemann problem 4 
from lZMOrj . WHAM is able to accurately resolve the contact dis- 
continuity and the height of the built-up state for this flow that 
has a large shearing velocity. 

cedure to decrease the order of the scheme to first order 
in sh ocks, does not produ ce post-shock oscillations for this 
test |Morsonv et al.|[2006l ). However, this flattening proce- 
dure requires parameter tuning for strong shocks (see sec- 
tion onij- 

5.17 ID relativistic Riemann problem 4 

The two initial states in this test have a large shearing ve- 
locity (table [TJ. Shearing flows are intrinsically hard in the 
relativistic case because of the coupling between the direc- 
tions through the Lorentz factor. This configuration, how- 
ever, does not lead to problems. Figure [TS] shows that the 
width of the contact discont inuity is the same as that of 
the WENO code from[ZM06; in our case, without the use 
of eigenvector decomposition, the built-up state near the 
contact discontinuity slightly undershoots compared to the 
analytic value. 




0.2 0.4 0.6 0.8 1 
x 

Figure 19. Density distribution for the one-dimensional relativis- 
tic Riemann problem 5 from lZMOrJ . For large transverse velocities 
as in this problem, one has to use an increased resolution in order 
to obtain the correct positions of discontinuities. 

discontinuities incorrectly. In order to get reasonable conver- 
gence to the analytic solution, one has to use an increased 
resolution for this test. Our code's re sult agrees with that 
of other codes at the same resolution (|Zhang fc MacFadvenl 
120061 : iMorsonv et al1l2006h . 

5.19 ID relativistic Riemann problem 6 

This is a very stringent test that probes the ability of a 
code to handle flows at extremely large Lorentz factors with 
strong shocks. It is a generalization of the nonrelativistic 
one-dimensional Noh problem (section 15. 5 \ in which a highly 
relativistic flow with 7 w 7 x 10 4 hits a reflecting boundary 
(see table [T}. Figure l20l shows that a reverse shock forms at 
the boundary x = 1 and moves to the left leaving the matter 
behind at rest. 

WHAM produces minimal post-shock oscillations com- 
parable to other schemes, such as a WENO code that uses 
eigenvector decomposition (ZM06). In the post-shock region 
the density reaches a maximum relative error of 1.5% in 
the vicinity of the re flecting wall. This is twice as small as 
with WENO (|ZM06l ) and HARM. Unlike the PPM algo- 
rithm l|ZM06h , we do not have to perform any fine-tuning of 
our code parameters for this test. 

A simila r yet more extreme test is discussed 
bv lAlov et all l|l999h for which they chose 7 = 2.24 x 10 5 , 
p = 7.63 x 10 -6 , and a resolution of 200 grid cells. As shown 
in Appendix IA10I their test problem can be barely resolved 
by a computer with a machine accuracy of ~ 10~ 16 (i.e. dou- 
ble precision). In order to minimize post-shock oscillations, 
they chose a Courant factor of 0.1 and tuned their recon- 
struction parameters. Using WHAM's standard numerical 
settings (e.g. Courant factor of 0.5), we obtain an error of 
< 0.5% in all quanti ties for this test, an error similar to that 
of lAlov et al.l (| 1999h (here we use their definition of relative 
error). The errors are dominated by the numerical solution 
having a shock width of a few points. 



5.18 ID relativistic Riemann problem 5 

This test is very similar to the previous one in terms of the 
initial conditions (table [1]). In this case, however, there is no 
shear; rather, both states are moving at a relativistic velocity 
in the transverse direction. Figure [19] shows that at a mod- 
erate resolution of 400 cells WHAM gets the positions of the 



5.20 2D relativistic shock-tube problem 

The initial conditions for this problem are shown in table [7] 
The test is run at a resolution of 400 x 400 until £f = 0.4. 
The contour plot of density is shown in figure 1211 This is a 
difficu lt highly relativistic two-dimen si onal Riemann prob- 
lem ( Del Zanna fc Bucciantinil 120021: iLucas-Serrano et alJ 
l2004l ; lzhang fe MacFadve n 2006; Morsonv et al.ll2006h . The 
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Figure 20. Density distribution for th e one-dimensional rela- 
tivistic Riemann problem 6 from ZM06. The numerical solution 
is shown with connected dots, and the analytical solution jZM06ft 
is shown with a solid line. WHAM is able to evolve extremely rel- 
ativistic supersonic flows with strong shocks with minimal Gibbs 
phenomenon. 



Table 7. Initial condi tions for the two-dimensional relativis- 
tic shock-tube problem l|Del Zanna fc Bucciantinil l2002'l. sec sec- 
tion ^. 20l The notation is the same as in table [4] The problem is 
computed on the domain (0, 1) X (0, 1) at a resolution of 400 X 400 
until time t F = 0.4 with T = 5/3. 
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Figure 21. 2D relati vistic Riemann prob- 
lem jDel Zanna fc Bucciantinil |2002| ; IZhang fc MacFadvenl 
120061) at the final time t F = 0.4, see section [5T201 The plot 
shows 30 equally spaced contours of log 10 p that go from 
-2.241 to 0.8243 in steps of 0.1057. WHAM is able to handle 
highly relativistic multidimensional flows with minimal Gibbs 
phenomenon. 



Lorentz factor reaches values larger than 25 inside the sharp 
diagonal jet-like feature in the lower left quadrant. The 
amount of oscillation in the lower -left quadrant is less than 
with the relativistic FLASH code (jMorsonv et al.1120061 ) and 
the RAM code (|ZM06t ), which uses finite-difference fifth or- 
der WENO with field-by-field decomposition. WHAM does 
not use such decomposition. The resolution of the diagonal 
feature in the upper-right quadrant is comparable to RAM. 

5.21 2D relativistic jet in cylindrical geometry 

We study a two-dimensional relativistic jet in cylindrical 
geometry. This problem couples many elements, which were 
separately tested in the individual problems described so far, 
in a real astrophysical setting: relativistic, highly supersonic 
flow, containing strong relativistic shocks, shear flows, and 
instabilities. Along with being of astrophysical significance, 
this test allows us to benchmark our numerical scheme 
against others. For the sake of comparison, we have u s ed the 
sam e conditions f or the test as IZhang fc MacFadvenl (|2006t ) 
and iMarti et al.l <|l997l , model C2). The size of the compu- 
tational domain is (0, 15) x (0, 45) in the (R, z) plane, and 
we use a resolution of 384 x 1152. The ambient medium is 
initially at zero velocity, with density p m — 1 and pressure 
Pm = 0.000170305. The jet is injected along the z-axis from 
the lower boundary with a density of pb = 0.01, pressure 
equal to the ambient pressure, pb = p m , and a speed of 
vt — 0.99c, which corresponds to a Lorentz factor of 7 ~ 7. 
Numerically, we implement the injection by assigning the 
values of the boundary grid cells within R < 1, z < to the 
state of the jet material. We use the zero-derivative outflow 
boundary conditions at the low-z (R ^ 1), high-z, and high- 



R boundaries and use the appropriate boundary condition 
on the axis R = 0. We evolve the system until tF = 100. 

The interaction of the jet material with the ambient 
medium forms an expanding bow shock and a Mach shock 
with a contact discontinuity in between that goes Kelvin- 
Helmholtz unstable, see figure 1221 The jet contains internal 
shocks, backflows, and shear flows. The jumps in the 7- factor 
on the right panel of the figure correspond to relativistic 
shocks crossing the axis of the jet. The essential structure of 
the jet, its head position, the shape of the bow shock, and the 
development of the Kelvin-Helmholtz instability agree with 
other codes' results (|Marti et al.ll 19971 ; IZhang fc M acFadven 
l2006| y 

5.22 Bondi flow in Schwarzschild geometry 

We study the order of convergence of our the scheme for sta- 
tionary sph erically symm etric accretion on to a non-spinning 
black hole l|Bondilll952r ). Even though the exact solution is 
spherically symmetric, for WHAM in 2D the Bondi prob- 
lem is in fact two-dimensional: in the f5-momentum equa- 
tion, a gradient of momentum flux (de(p g sin 9) in the non- 
relativistic limit) is balanced by a source term (p g cos 9) 
that numerically cancel out to within the truncation er- 
ror ijGammie et al.l 120038 ) . We find that in order to obtain 
an accurate solution, one has to properly average the source 
terms on the r.h.s. of the equations of motion ([6]). Fur- 
ther, for this test WHAM uses Kerr-Schild coordinates that 
have nonzero space-time mixing even for the case of a non- 
spinning black hole: 7 out of 10 components of the metric 
are nonzero, so this problem involves many of the terms 
that appear in the general equations of motion. For this test 
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Figure 22. Rclativistic two-dimensional jet problem in cylindri- 
cal geometry: the left panel shows the color-coded distribution of 
the logarithm of the fluid frame rest-mass density (red denotes 
high values and blue low values), and the right panel shows the 
^-dependence of 7 along the axis of the jet. T he sim ulation was 
run at a resolution of 384 X 1152, the same as in ZM06, and the re- 
sult is mirrored across the axis in order to obtain the figure in the 
left panel. For more detail, see section [5.211 WHAM reproduces 
all essential elements of the jet structure for this astrophysically 
relevant problem. 



problem we set the source term analytically in order to test 
convergence near machine precision accuracy. 

We follow t he se t up of the pro b lem as discussed 
by lHawlev et al.l (| 19841 ); iGammie et all (|2003l ). We fix the 
radius of the sonic surface, R s — 8, the adiabatic index of the 
gas, r = 4/3, the gas adiabat, K = 1 (= p g /p F ), and choose 
the mass of the black hole, M = 1. These determine the 
critical solution and the three integrals of the problem: the 
radial flux of mass (= y/—gpu r ), the flux of energy (= T r t ), 
and the entropy of the flow. The solution can be found semi- 
analytically by solving an 8t h order polynomial e quation for 
the temperature of the flow l|Hawlev et al.lll984l ). 

We treat the boundary conditions as follows. We use 
the value of the 'baryon flux' Ci from the outermost active 
grid cell together with the critical values of K and the 'en- 
ergy flux' C2 to set the state of the outer-r boundary cells 
(for definitions of the integrals G\ and C2, see lHawlev et al] 
1 19841 ). At the inner-r boundary we choose the states of the 
cells to correspond to the three integrals of the flow deter- 
mined from the numerical solution at the innermost active 
grid cell in the active grid. At the ^-boundaries we use the 
usual polar axis boundary conditions. As the system ap- 
proaches steady state, the mass flux is determined by the 
critical condition at the sonic surface. Since in steady state 
the fluxes become uniform in space within the truncation 
error, extrapolating the flux allows us to approximate the 



Table 8. Relative L\ -errors in density and the order of conver- 
gence of the WHAM, HARM, and WHAM-IS schemes for the ID 
and 2D spherical accretion problems (section 15. 22 H , at TV X N res- 
olution. At the same resolution, the higher order WHAM scheme 
provides a much smaller error than the other two schemes. The 
averaging of the source terms is crucial for preserving the high 
order accuracy. Radial velocity and internal energy similarly con- 
verge to the analytic solution. 



N 


16 


32 


64 


128 


256 


512 


WHAM 
Order 


1.6e-5 


ID Bondi problem 
1.8e-7 1.5e-9 1.6e-ll 
6.4 6.9 6.6 


1.2e-12 
3.7 


5.0e-15 
7.9 


TT A Tl TV iT 

HARM 
Order 


9.7e-4 


2.1e-4 
2.2 


4.9e-5 
2.1 


1.2e-5 
2.0 


2.8e-6 
2.0 


6.9e-7 
2.0 


WHAM-IS 
Order 


2.0e-3 


5.2e-4 
2.0 


1.3e-4 
2.0 


3.2e-5 
2.0 


7.9e-6 
2.0 


2.0e-6 
2.0 


WHAM 
Order 


1.6e-5 


2D Bondi problem 
1.8e-7 1.5e-9 1.6e-ll 
6.4 7.0 6.5 


1.2e-12 
3.7 


4.9e-15 
7.9 


HARM 
Order 


6.4e-4 


1.6e-4 
2.0 


3.9e-5 
2.0 


9.7e-6 
2.0 


2.4e-6 
2.0 


6.0e-7 
2.0 


WHAM-IS 
Order 


1.4e-3 


3.3e-4 
2.0 


8.2e-5 
2.0 


2.0e-5 
2.0 


5.1e-6 
2.0 


1.3e-6 
2.0 



boundary conditions very accurately, compared to extrapo- 
lating, e.g., the nontrivially varying density. This also avoids 
the order of the extrapolation limiting the order of conver- 
gence of the scheme. 

We have performed the problem in ID and 2D in modi- 
fie d Kerr-Schild c o ordin ates, the same coordinate system as 
m bammic et al] (|2003l ). over the domain r £ (1.9,20) in 
units of GM/c 2 . We initialize the problem with the analytic 
solution and let the system evolve for t = 1000 in units of 
GM/c 3 , corresponding to several sound crossing times, by 
which time the system reaches steady state. We then com- 
pute the relative Li-error norm in density between the final 
solution and the initial conditions over the active grid. The 
results are shown in table [8]for various resolutions. WHAM 
converges at an order higher than 5 and makes a much 
smaller error than the WHAM-IS scheme, which does not 
average the source terms, or the HARM scheme. The lat- 
ter two schemes converge at second order as expected, since 
they do not properly account for the difference between the 
point and average source terms and/or fluxes and conserved 
quantities. The nonuniform convergence rate appears to be 
due to the sensitivity of the monotonicity indicator (see sec- 
tion \KE§ for particular parts of the solution. 

5.23 Equilibrium torus in Kerr geometry 

Finally, we consider a rotating fluid torus surrounding a 
spinning black hole that is in equili brium under the ac- 
tion o f pressure and centrifugal forces (jFishbone fc Moncriej 
Il97q ). We choose the black hole spin to be a/M = 0.95, 
and take the following parameters for the torus: loca- 
tion of the inner edge rj n = 3.7, angular momentum per 
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Table 9. Relative Li-error in density and order of convergence 
in density of the WHAM and HARM schemes for an equilibrium 
torus problem, section [523] at N X N resolution. Similar conver- 
gence is observed in all other quantities. For this general relativis- 
tic problem in the Kerr metric, WHAM achieves asymptotic fifth 
order convergence. 



N 


16 


32 


64 


128 


256 


512 


WHAM 


1.6e-2 


1.7e-3 


2.3e-4 


2.8e-5 


1.4e-6 


2.7e-8 


Order 




3.3 


2.9 


3.1 


4.3 


5.7 


HARM 


1.5e-2 


3.5e-3 


6.6e-4 


1.4e-4 


3.3e-5 


8.3e-6 


Order 




2.1 


2.4 


2.2 


2.1 


2.0 



unit mass tt'wy = 3 .85 (as in iMcKinnev fc Gammi el 120041 : 
iGammie et al.l 2004), and equation of state p g = Kp F with 

K = 10" 3 , r = 4/3. 

The exact iFishb one fc Moncriej torus solution is em- 
bedded into a vacuum. In order to avoid having zero 
densities that are problematic numerically, we embed 
the torus in a dynamically unimportant atmosphere by 
introducing floors on the density and internal energy: 



-3/2 



= 10- 6 (r/r fa )- 5 / 2 (as 



Pmin = 10~ jr/Hp 

IGammie et"al1l2003h . 

We use th e same modified Kerr-Schild coordinates as in 
IGammie et~ai] (|200St ). concentrating the resolution toward 
the midplane, and run the simulation until tp = 10 in the 
2D domain (r,6) G (0.98r h , 20) x (— k/2, 7r/2), where r h 
is the radius of black hole horizon. We compute the relative 
L\ -error norm between the final solution and the initial con- 
ditions over the region where p > 0.02/9 max to minimize the 
influence of the atmosphere on the convergence results. The 
results are shown in table [5] We see that WHAM asymptot- 
ically converges at fifth order0 while HARM converges at 
second order. 

This final example tests all the general relativistic as- 
pects of WHAM - equations of motion, metric, connection 
coefficients, and source terms - in Kerr space-time. 



6 LIMITATIONS 

The hydrodynamical numerical scheme we have described 
actually does include magnetic fields, with the divergence- 
free constraint ()5d|) applied using the FL UX-CT constrained 
transport method jGammie et al.l l2003). This method aver- 
ages the fluxes of the magnetic field in such a way that 
the associated update to the magnetic field is guaranteed to 
preserve a certain numerical representation of magnetic field 
divergence. However, the magnetic fields are treated to lower 



8 At low resolutions, the truncation error of WHAM is dominated 
by unresolved edges of the torus where the conserved quantities 
rapidly fall off to zero and the reconstruction becomes second 
order. Even though the torus is resolved with about 30 grid cells 
at a resolution of 64 X 64, the drop of Uq = yj—gpu 1 at the 
torus edge is resolved with only 4 grid cells. Because of this, at 
low resolutions HARM is comparable to WHAM, but at high 
resolutions WHAM is far more accurate. 



order than the hydrodynamical quantities. In a followup pa- 
per (McKinney, Tchekhovskoy, & Narayan in prep.), we will 
describe a consistent high order scheme that generalizes our 
hydrodynamical method to a full MHD method with a con- 
strained transport method that is a high order version of the 
stagg ered grid method used by Athena l|Gardiner fc Stonel 
120051 ). We expect that the stiff regime associated with the 
MHD equations in the highly magnetized limit will signifi- 
cantly benefit from our higher-order reconstruction method 
since the effective magnetic Mach number is M mag , which is 
on order of 10 or larger for models of black hole or neutron 
star systems. 

WHAM is capable of performing simulations in three di- 
mensions. However, a description of our method and a suite 
of 3D tests will be given in a followup paper (McKinney, 
Tchekhovskoy, & Narayan in prep.). 

We have not implemented any performance optimiza- 
tions in our scheme since much has been already optimized 
in HARM which WHAM is based upon. Currently, for one- 
dimensional problems, WHAM i s about 4 t i mes slo wer than 
the third order method used bv lMcKinnevI (|2006bl ). We ex- 
pect that further optimization will lead to a factor of sev- 
eral speed increase and make WHAM competitive with any 
scheme of similar accuracy/order. Our experience suggests, 
however, that even without optimizations the duration of 
simulations will not be the bottleneck in completing science 
projects. 



7 PROPOSED APPLICATIONS 

The numerical scheme that we have developed provides su- 
perior accuracy in highly supersonic flows as well as flows 
in which the energy scales are very different. The ap- 
plications where the scheme is particularly advantageous 
therefore include the study of supersonic disk winds, jets, 
and pulsar magnetospheres. There have bee n numerical 
studies of t hese systems within the force-free ([Komissarovl 
20021. 12006J ; lMcKinnevll2006al; IMcKinnev fc Naravanll2007 



McKinnevf 12006c ; ISpitkovskyl 120061 ; iNaravan et all l20o' 
and full MHD (jKomissaroyl I2005T IMcKinnev fc Gammiel 
20041; iBucciantini et al.l l2006; Mc Kinney fc N arayan 2007al : 



McKinnevH2006bl ) limits. IMcKinnev! (|2006bl ) achieved high 
magnetization, b 2 / (p + u g +p g ) ~ 10 3 , and could handle val- 
ues of up to 10 4 but it was not clear whether the evolution 
was accurate in that regime. However, even such high mag- 
netization values are still below the actual observed values 
of 10 5 - 10 6 in the case of pulsar outflows. 

We plan to study jet and pulsar systems in the limit of 
cold GRMHD, i.e. with the internal energy of the plasma 
vanishing exactly. This limit should allow us to achieve 
higher magnetisation factors more easily. Since our code 
does not require the use of eigenvector decomposition, only 
the equations of motion have to be changed to study systems 
in this limit. 

Our code can also be used to perform three-dimensional 
studies of thin discs around black holes that involve highly 
supersonic motion in the ^-direction. Such studies would 
determine how jet speed and power is related to disk 
thickness, address the problem of the multidimensional 
stru cture of disks near the innerm ost stable circular or- 
bit l|Beskin fc Tchekhovskovl l2005h , and test the applica- 
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bility of well-known analytic models (|Shakura fe Sunvaevl 
11973 : iNovikov fc Thornelll973h there. 



8 CONCLUSIONS 

We have developed a high order accurate conservative finite 
volume general relativistic hydrodynamic (GRHD) scheme 
called WHAM that is robust in the vicinity of shocks 
and is very accurate in smooth flows even for an arbi- 
trary metric and coordinates. It compares very well with 
its special relativistic hydrodynamic analogs and is, as 
far as we know, the first GRHD code that converges at 
fifth order in smooth flows. In con t rast to most high or- 
der schemes ll Jiang fc Tadmorlll998l: iTitarev fc Toroll20o3 : 



iFene : et al.l [20041 : IZhang fc MacFadvenl l2006h . WHAM uses 
primitive quantities to reconstruct the values of quantities 
at interfaces. This is a safe choice in the relativistic regime 
because it guarantees that the interface values, which deter- 
mine inter-cell fluxes, are physical. 

WHAM avoids the use of the computationally intensive 
eigenvector d ecomposition approach used by many high or- 
der schemes ll Jiang fc Tadmoij [l998l: ITitarev fc Tordl2004l; 
Feng et all 12004 IZhang & MacFadvcn 200l IXing fc Shul 



2006). Instead, we adaptively reduce the order of the scheme 



near discontinuities, while, unlike the standard WENO for- 
malism, avoid excessive reduction in smooth monotonic 
flows. The resolution of contact discontinuities we obtain is 
comparable to that of WENO schemes that use eigenvector 
decomposition without contact steepeners. 

Unlike finite differ e nce schemes (e.g., Jiang fc Tadmorl 
1 19981 : iFeng et all 120041 : IZhang fc MacFadvenl 120061 ). which 
conserve the integrals of motion up to truncation error, our 
finite volume scheme conserves them exactly, to machine 
precision. Our scheme performs the proper conversion be- 
tween the average and point values of conserved quantities, 
fluxes, and source terms using fifth order WENO-type recon- 
struction. We find that the alternative me thod of quadra- 
tures is less accurate for this purpose (e.g ., ITitarev fc Torol 
120041 : iNoelle et alj|2006l : IXing fc Shull2006l ). 

We have demonstrated that our scheme provides accu- 
rate solutions in highly supersonic flows, which are intrin- 
sically hard for conservative numerical schemes. In particu- 
lar, the scheme is able to accurately evolve shocks for which 
both pre-shock and post-shock material moves supersoni- 
cally through the computational grid. The scheme is likely 
to be particularly good for studies of relativistic and/or su- 
personic flows near compact objects. In the future we plan 
to extend the formalism to include magnetic fields and study 
models for which the magnetization is large to identify why 
present schemes have difficulties in such regimes (McKinney, 
Tchekhovskoy, & Narayan in prep.). 
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APPENDIX A: IMPLEMENTATION NOTES 

Al Truncation error analysis for a smooth high 
Mach number flow 

This section presents a detailed error analysis of the evo- 
lution of the Hubble-type diverging flow, see sections 13.21 
and l5.ll We are particularly interested in the consequences 
of neglecting the difference between the cell average and cell 
centred values of conserved quantities. As we show below, 
this leads to a secular error in the evolution of the internal 
energy so that it becomes negative. The error is first order 
in time and second order in space, independent of the actual 
order of discretisation of the scheme in space and time. This 
means that lowering the time step does not diminish the 
error; only by significantly increasing the spatial resolution 
can the error be reduced. 

Firstly, we verify that solution (|12[) satisfies the conser- 
vation laws: 



dp 
dt 

djpv) 
dt 

dUg_ 
dt 



v po 



{i + v' ty 
(i + v> ty ~ 



9(pv) 
dx ' 

dx 

_ d[(u g +p 3 )v] 
dx 



(Al) 
(A2) 
(A3) 



so the system (|12[) correctly describes the time evolution 
of the initial distribution described by setting t = in the 
system. 

Let us assume that the spatial reconstruction on the 
primitive quantities is exact but one does not take into ac- 
count the difference between cell averaged and cell centred 
values of conserved quantities. We use the Euler method to 
discretise time stepping and check if the method converges 
at second order in time as one would naively expect. We only 
need to consider the first time step since for all successive 
time steps the result of the previous time step(s) can be con- 
sidered an initial condition. We define a uniform numerical 
grid consisting of grid cells Ai — [a^i-i/2, ^i+1/2] of size h. 
For spatial discretisation we use the simplest finite volume 
scheme: we linearly interpolate the values of primitive quan- 
tities from cell centres to cell interfaces within each of the 
grid cells. Since in the above solution all primitive variables 
depend linearly on x, this linear interpolation is exact. 

Consider the ith grid cell, Ai. For clarity we shall omit 
the index i where it does not lead to confusion. Let us in- 
tegrate the equations from i = to t = At. The conserved 
variables are the density p, specific momentum p = pv, and 
the total specific energy E = pv 2 /2 + u g . The discretisation 
of conservation laws then takes the form: 

A<p) = - At [p v' (x + h/2) - p v' (x - h/2)] /h 

= - Po v' At, (A4) 

A(p> = - Ai [p v' 2 {x + h/2) 2 - Po v' 2 (x - h/2) 2 ] /h 

= - 2p v Q v' At, (A5) 
A{E) =-At [p v' 3 (x + h/2) 3 /2 - p v' 3 (x - h/2) 3 /2] /h 

~ At(Ug +Pg) v' 

= - (pv 2 /2) 3v' At (1 + 7i2 h 2 /x 2 ) - (u g +p g ) v' At, 

(A6) 

where angle brackets indicate averaging over the grid cell 
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and index '0' indicates values taken at t = 0. The latter 
equation corresponds to the following conservation law: 

(A7) 

The above expressions ()A4() - (|A6|I are so far exact up to first 
order in time. Now suppose we remove the angle brackets 
on the l.h.s. in the above expressions, i.e. do not make a 
distinction between cell centred and cell averaged quantities. 
Equations (|A4[l and (|A5I) will still be exact up to first order 
in time for the particular problem considered, i.e. they will 
still give the updates to cell centred quantities p and p that 
are first order accurate in time. This is because p and p — pv 
axe linear functions of x and thus their average value in any 
grid cell is the same as the point value at the cell centre. This 
is not the case, however, for equation ()A6|) where energy E 
is a nonlinear function of x: E 7^ (E). This also can be seen 
from the expansion of equation (|12|l , according to which the 
cell centred conserved quantities evolve according to 



Ap = p (-v' At + v' 2 At 2 -...), 
Ap = po (-2v' At + Sv' 2 At 2 -...), 
AE = (pv 2 /2) (-3v' At + 6v' 2 At 2 - 



(A8) 



) 



+ (u g + Pg )o (-v' At + (r + l)«o 2 At 2 /2 -...). 

So, by treating the update to the average of the energy 
in the grid cell as the update to the point value of energy, 
we are making an absolute error in AE equal to 



1 2 / L 2, h 2 



pov'oh 2 At. 



£(AE) = A{E)-AE = --povfohTAt-^ _ s< 

(A9) 

Note that up to first order in time, £(Ap) — £(Ap) — 0. 
Since the internal energy is computed from the values of 
conserved quantities via 



Ug — E 



(A10) 



(All) 



-p7(2p), 

the error in the update to internal energy is 
£(Au g ) = £{AE) = -\p v'oh 2 v' At. 

o 

It is first order in time for every step, which means that 
using smaller steps (i.e. lowering the Courant factor) does 
not lead to convergence. We note that using higher order 
discretisation in time and/or space does not help to avoid 
this error. If we were to use a higher order time stepping 
discretisation scheme, the time stepping error would remain 
the same because higher order time stepping schemes as- 
sume that each simple first order time step is accurate to 
first order. However, here this is not the case. So for higher 
order time stepping error (|Alip remains the same. Only by 
increasing the spatial resolution can one get convergence of 
second order in space, independent of the actual spatial and 
temporal order of approximation of the scheme. 

We can estimate the error in the internal energy at a 
characteristic time of the evolution by setting At — l/v' , 



£(Au g )\ 



p v' 2 h 2 



-Mi 



(A12) 



Uq \v' Q At=l Uo 

This error is on the order of the minimum of the Mach num- 
ber on the grid squared. Therefore, for large Mach numbers, 
the schemes that do not take into account the difference be- 
tween the point and average values of conserved quantities 



Table Al. Optimal weights for various types of W ENO-5 re- 
construction. According the convention of Shu (1997) adopted in 
this paper, do, di, and di are the right, central, and left optimal 
weights. 



Reconstruction type 


d 


di 


d 2 


Centre to left face 


1/16 


5/8 


5/16 


Centre to right face 


5/16 


5/8 


1/16 


Average to centre 


-9/80 


49/40 


-9/80 


Centre to average 


-17/240 


137/120 


-17/240 



make very large errors in the evolution of the internal en- 
ergy and require the use of the resolution proportional to 
the Mach number of the flow in order to correctly capture 
the evolution of the internal energy. For instance, for a Mach 
number M m i n = 100, the resolution has to be increased by 
about two orders of magnitude. See section I5TT1 for a numer- 
ical verification of these results. 



Reconstruction matrices and optimal weights 
for WENO-5 



A2 



In WENO-5 the integer k in section 14.51 is equal to three. 
Thus we have 3 stencils, each of length 3. The reconstruction 
matrices in equation (|18p for average to centre and centre 
to average reconstructions, c^ c ' and c^ a ' , are the inverse of 
each other: 



„( ac ) 



„( ca ) 




The reconstruction matrices for centre to left interface and 
centre to right interface reconstructions, c+P and cf- , are: 



(cl) 
c ■ ■ — 



15/ 

/s 


5/ 


7 


% 


A 


-7 


-7s 


3/ 

A 


7 


% 


3/ 

A 


-7 


-7« 


3/ 

A 


7 


7s 


5/ 

- A 





The optimal weights for centre to edge, average to cen- 
tre, and centre to average reconstructions can be computed 
by equating the expressions (|19[1 and (|20[) . For completeness, 
we give them here, see table IA1I Note that for average to 
centre and centre to average reconstructions, negative op- 
timal weights appear that require a special treatment. See 
section |4"7H for a discussion. 



A3 Choosing the value of e 



The standard WENO weights computation procedure uses a 
fixed val ue of e for the purpose of avoi ding division by zero 
23J) jjiang fc Shul 11999 : Isiiulll997r i. However, setting e 



© 2006 RAS, MNRAS 000. rflMl 



26 A. Tchekhovskoy, J. C. McKinney, & R. Narayan 



to a constant value independent of the function being inter- 
polated sets artificial, irrelevant scales for a problem. Any 
discontinuity in the interpolant, such that smoothness indi- 
cators for the stencils that cross that discontinuity are much 
smaller than e, is mistakenly treated by the reconstruction 
as part of smooth flow. 

Let us consider e <C 1. Then a contact discontinuity, 
which is a density jump with Ap = 1 at x = 0, 



p(x) 



2, x > 
1, x < 0, 



(A13) 



is correctly treated by the reconstruction, i.e. the reconstruc- 
tion avoids using the stencils that pass through the disconti- 
nuity. However, if the same density jump is recast in different 
density units, 



p(x) = 



2e, x > 
e, x < 0, 



(A14) 



it is not avoided by the reconstruction because in (|23[) 
e > /j r ~ (A/5) 2 ~ e 2 . The addition of e to the smoothness 
indicators in this case effectively hides the jump from the 
reconstruction and leads to an oscillatory reconstruction. 

We propose the following algorithm for the adaptive 
choice of e. We leave the weights computation (|23[1 - (|24p 
the same but modify the smoothness indicators that come 
into it: 



P' r = Pr +4v 2 || +5, 



(A15) 



where 



is the sum of v 2 within the WENO stencil and S 



is the smallest positive number the floating-point type used 
in the numerical implementation of the method can hold; it 
is added to avoid division by zero when ||w 2 || vanishes. We 
choose an e such that it is larger than the relative machine 
precision in the computation of smoothness indicators and 
at the same time is negligible compared to them in shocks, 
e = (ce machinc ) 2 . Here e mac hino is of order of relative machine 
precision for the float type used and c is a sufficiently large 
constant. We choose c = 100, so for double precision we have 

^machine ~ 10~ 15 ande=10- 26 . 

Further, in order to extend the dynamic range of the re- 
construction and avoid division by small numbers (smooth- 
ness indicators become small for small ||fJ 2 ||), prior to plug- 
ging p! r 's into the weights computation (|23[1 we rescale them 
so that they add up to unity: 



ft 



(A16) 



where ||#||=£^ -1 #. 

Finally, we use the usual weights computation proce- 
dure (|23[) - (|24p in which we use modified smoothness indi- 
cators /?" (|A16[) and the value of e determined above. 



A4 At what order does WENO-5 actually 
converge in smooth flows? 

Assuming we can neglect e in the definitions of the weights, 
then in order to comply with (|21[) we need to have for 
r = 0, k- 1: 



p r = D [l + 0(/ l fc - 1 )] 



where D is a positive constant common to all stencils at a 
given resolution (so D may depend on the grid cell size, h). 

Let us verify if for a smooth function v(x) smoothness 
indicators (|22[) satisfy requirement (]A17[) for the case of 
WENO-5, i.e. k — 3. The smoothness indicators for cell 
centre to interface, average to centre, and centre to average 
reconstructions are the same for the case of WENO-5 (see 
appendix I A6|) and are given by: 

1 2 13 2 

Po = - (~3vi + 4v i+1 + v i+2 ) + — (vi - 2v i+ i + Vi+2) , 



4 
1 
4 
1 

Using Taylor expansion, this can be cast into 

#> = \{2v'h - 2/3v"'h 3 - l/4v W h 4 ) 2 + ^(v"h 2 + v"'h 3 ) 2 , 



Pi = 7(^-2 — 4ui-i + 3vi) 2 + ^-(vi-2 — 2«»_i + Vi) 2 . 



(A19) 



( 3 1 = ± { 2v'h+l/3v"'h 3 ) 2 + 1 ^(v"h 2 ) 2 , 

p 2 = h 2v 'h - 2/3v"'h 3 + l/4v (4) h 4 ) 2 + ^(v"h 2 - v'"h s ) 2 , 

with the truncation error 0(h 6 ). Here the derivatives are 
evaluated at the cell centre x = x*. i/™' = v'™' (ajj). If v' 7^ 0, 
then 

f3 r = {v'hf [l+0(h 2 )] , r = 0,l,2, (A20) 

which satisfies (|A17|) and agrees with Ijiang fc Shul (| 1996h . 
This means that if the ^'-containing terms dominate in 
(|A19I) . then the smoothness indicators give nearly optimal 
weights that satisfy (|2ip. and hence the truncation order of 
WENO reconstruction is 0{h 2k - 1 ) ~ 0(h 5 ) . 

However, at simple extrema with v' = and v" 7^ 0, or 
more precisely with \v'\h <C |v'"|/i 3 \v"\h 2 , we have for 
smoothness indicators, 



Pt = ^iv"h 2 ) 2 [1 + 2 (v"'/v")h + 0(h 2 )] , r = 0, 1, 2, 

(A21) 



(A17) 



which means that neither (|A17[l nor (|21|> are satisfied, and 
WENO reconstruction (I19|) has a larger truncation error, 
(D(h 2k ~ 2 ) ~ 0(h 4 ), near the such extrema. This reduction 
in the order of the scheme occurs because the second deriva- 
tive approximation in the smoothness indicators (|A18[) is 
only first order acc urate. This disagrees with the result of 
Ijiang fc Shul |l996). They erroneously consider the numer- 
ical approximation to the 2nd derivative to be 2nd order 
accurate, i.e. they are missing the term proportional to v'" 
in (|A21[) . Therefore in general the approximation order of 
the WENO-5 scheme at extrema reduces to 4, and only in 
the special case of extrema with v'" = does WENO-5 give 
fifth order approximation, e.g. for v(x) = sin a;. 

Finally, if both first and second derivatives become 
small compared to higher order ones (at a high order ex- 
tremum, a critical point, or in a smooth region of the 
flow dominated by high-order derivatives) the order of the 
WENO-5 scheme reduces to the third one. This happens 
in contact discontinuities and leads to their excessive dif- 
fusion, this also happens near inflection points, e.g. of the 
v(x) — x s function: here WENO-5 is unable to capture the 
dependence near x — exactly even though a fifth order 
polynomial would be able to. 
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A5 High order approximation in smooth 
monotonic flows 

WENO-5 scheme aims to obtain the flattest reconstruction 
profile inside a grid cell of interest: it favors the stencils that 
minimize the absolute value of the first and second deriva- 
tives within the grid cell. In smooth flows that are locally 
well approximated by a parabola these derivatives do not 
change significantly between the stencils, so such a prescrip- 
tion leads to nearly optimal weights given to stencils and 
results in the fifth order of approximation. However, if any 
terms higher than second order ones are significant in the 
Taylor expansion of the function, the WENO-5 scheme re- 
duces the reconstruction order down to the third one. This 
reduction leads to unwanted diffusion of contact discontinu- 
ities that becomes very large without the use of field-by-field 
decomposition. 

To avoid such a reduction of the order of the scheme, 
we use a fifth order interpolation polynomial for reconstruc- 
tion if the polynomial and its derivatives are monotonic. 
This minimizes the amount oscillation in the reconstructed 
solution and maintains a high order of interpolation in the 
monotonic regions of the flow. We use this approach for the 
centre to edge reconstruction. A similar idea is used in the 
PPM reconstruction procedure (|Colella fc Woodward|[l984l ; 
iMignone et al.ll2005h . 

Consider an exponential dependence f(x) = e x with 
the grid cell centres located at ii = hi, h = 1. In 
this case the WENO-5 weights (|24|l are (ujo, uii,U2) oc 
(0.007do, 0.09di , 0.9^2) independent of the value of i; here 
dr's are the optimal weights (see text after equation (1201) ). 
WENO-5 is giving about 90% of the weight to one of the 
stencils, therefore the WENO-5 reconstruction order is ef- 
fectively reduced from 5 to 3. In this case, the fifth order 
polynomial interpolation through, e.g., points X-2, ...,X2 
would have no oscillations at the interval [x-i, xi] and could 
well be used. 

Therefore, prior to applying the standard WENO-5 re- 
construction inside an ith grid cell, we consider the fifth or- 
der reconstruction polynomial fs{x). It interpolates the set 
of discrete values fi-2, ■ • ■ , fi+2 at points Xi-2, . . . , x i+ 2- We 
check if all of the polynomial's derivatives are monotonic at 
the interval [xi-\, Xi+i] and if they are, we use this polyno- 
mial for interpolation instead of the WENO-5 interpolation. 
For a fifth order polynomial, it suffices to check that each 
of the first three derivatives, f$"\ n = 1,2, 3, has the same 
sign at both ends of the interval; below for compactness we 
shall omit the index '5'. 

We now describe an indicator that vanishes unless the 
function values fi-2, ■ ■ ■ , fi+2 form a monotonic sequence 
and the interpolation polynomial fix) and all of its deriva- 
tives are monotonic at (xi-i, Xi+i) (in which case it is 
unity). In order to avoid switching between the schemes, 
we design the indicator to be a smooth function of dis- 
crete function values. We define the minimum value by 
absolute magnitude of an nth derivative at the interval, 

)j . In order to en- 



/£& = MINMOD [/< n >(zi-i),/ {n) 0= 
sure monotonicity of discrete values, we redefine f t 
MINMOD [f£l, - fi-2, fi+2 - fi 
pute the minimum absolute 



(i) 



I . Finally, we corn- 
value of nondimensionalized 



min \f^ilh n . We refer to the weight that we give to the 

n=l, 2,3 1 

fifth order polynomial reconstruction as the monotonicity 
indicator, 



:{o,mm[l,f^J(V~e\\f\\+S)]}, (A22) 



where ||/|| = X^ii^l/il * s the norm °f f( x ) a t t ne stencil 
and 8 - added to avoid division by zero - is the minimum 
positive value that a floating-point variable can hold on our 
machine. In the above formula e controls how fast the use of 
the fifth order polynomial for reconstruction kicks in when 
the function and its derivatives become monotonic. For con- 
sistency, we choose e to be the same as the one in the WENO 
weights computation procedure (|A15I) . (|A16[) . (f2"3")l . P4")) . 

The reconstructed value we use is given by a linear com- 
bination of the fifth order accurate value of the polynomial, 
fs, and that of the WENO- type reconstruction, /weno: 



/ = q/ 5 + (1 — a)/wBNO, 



(A23) 



where we have returned back to the notation /s for denoting 
the fifth order reconstruction polynomial. 

Since the reconstruction due to the fifth order poly- 
nomial is equivalent to the WENO-5 reconstruction with 
weights set to optimal ones, equation (|A23|I is equivalent 
to using the WENO-5 reconstruction with modified weights 



io' r = ad r + (1 — a)co r , 



,k - 1. 



(A24) 



This provides a smooth transition between using the poly- 
nomial and the WENO-type reconstructions. 



derivatives of the reconstruction polynomial f(x), f 



CO 



A6 Higher order smoothness indicators 

High order WENO-type schemes provide better convergence 
properties near critical points of smooth flows than their 
lower -order counterpar ts but have a higher computational 
cost |Qiu fc !§h3[2002). In this section we summarize the 
properties of higher order smoothness indicators. In particu- 
lar, we concentrate on the convergence properties of WENO- 
type schemes near the critical points of smooth flows: when 
does a WENO-n scheme provide the maximal, nth, order of 
convergence? 

One can show for any k (at least up to k = 7 as we have 
checked in Mathematica) that any smoothness indicator (|22|l 
can be cast in the form analogous to that of (|A19|) : as a linear 
combination of the squares of discrete approximations to the 
derivatives at x = Xi. For instance for the case of k = 4, we 
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have for the smoothness indicators: 

00 = — 732 (— 3vi + 4vi+i — v i+2 ) 2 

+ 7s2(— + 18«»+i - 9u l+2 + 2-y i+3 ) 2 

+ 13 /i2(2ui — 5vi+i + 4v i+2 — fi+3) 2 

+ 316 %88o(— Vi + 3v i+ i — 3v i+ 2 + v i+a ) 2 , (A25) 

01 = + l /l${Vi+l — Vi-i) 2 

+ 1 /4s(—2vi-i — 3vi + 6v i+ i — Ui+2) 2 

+ 13 /l2{vi~l — 2Vi + v i+ i) 2 

+ 3109 /288o(— Vi-i + 3vi — 3vi+i + v i+2 ) 2 , (A26) 

fh = + 7l6(«»+l — «i-l) 2 

+ 7-4s(Wi-2 — 6«j_l + 3«; + 2w i+ i) 2 
+ 13 /i2(w»-i - 2Vi + v i+ i) 2 

+ 3109 /2880(-Ui-2 + 3ttj_l - 3Wi + v l+1 ) 2 , (A27) 
03 =~ l /z2(Vi-2 - 4Vj_l + 3fi) 2 

+ 732(-2Vi-3 + 9Ui_2 - 18Ui_l + ll«i) 2 

+ 13 /i 2 (-^-3 + 4^-2 -5«i-i +2u t ) 2 

+ 316 72880(-^-3 + 3l?i-2 - 3«i_l + Vif. (A28) 

This is a much more conceivable representation of 
higher order smooth ness indicators than the one used in 
iBalsara fc Shul l|2000l ). In fact, this form allows one to easily 
prove the properties of higher order smoothness indicators 
the same way as we did for the case k = 3. Again, using 
Taylor expansion, we get: 

O = -- (v'h-h^h 3 - I„«A 4 Y + H (v'h+ i,Wy 
8V 3 4 J 8\ 4 J 

+ H (V'A 2 _ l± M h A\ 2 3169 / (3) a 3 v (4) h 4y 
12 V 12 / 2880 V 2 / ' 

(A29) 

fh= + \ (v'h + i^ 3 ^ 3 ) 2 + \ fhv' - -^V 4 >) 2 

+ H (V'ft 2 + L v W h A 2 + 31M / (3) a I ft 4 ( 4)\ 2 
12 V 12 / 2280 V 2 7 ' 

(A30) 

02=- I [v'h + JtW) 2 + \ (t/h + ^<*V) 2 

+ H ( w v + L v w h A\ 3109 / (3) 3 _ i ft4 (4) y 

12 V 12 / 2880 V 2 / ' 

(A31) 

3 = - i (V A - i« (3 V + i« (4) ft 4 ) 2 + I (»'A - \hWj 
+ i 3 . (V'A 2 _ H^(4)\ 2 3169 / (3) 3 _ 3 4 ( 4) \ 2 

12 v 12 y 2880 v 2 y ' 

(A32) 

with the truncation error of 0(v'h 6 ) + 0(v"h 7 ) + 0(h s ). 
If v' 7^ 0, we get the analog of equation (|A20|I . 

r = {v'h) 2 [1 + 0(h 2 )] , r = 0,l,2,3; (A33) 

further, for v = 0, v" / Owe obtain the analog of equation 

(EE), 

/3 r = 13/12(u"/i 2 ) 2 [1 + £>(h 2 )] , r = 0,1,2,3. (A34) 



This shows that the smoothness indicators for WENO-7 
(with k = 4) are capable of handling the extrema, where 
v = 0, without any reduction in order: the order of the 
scheme is then 0{h' 2k ~ 1 ) ~ 0(h 7 ). However, when both 
v' = and u" = 0, 

/3 r = 1043/960(t/ 3) /i 3 ) 2 [l + C(A)], r = 0,1,2,3, (A35) 

so the order of the scheme reduces to 0(h 2k ~ 2 ) ~ 0(A 6 ) for 
this case. 

More generally, consider a WENO-(2fc — 1) scheme, 
which uses stencils of length k. As can be verified for each 
particular k, the numerical approximations to v' , . . . , i/* -2 ' , 
that the smoothness indicators are comprised of, have trun- 
cation errors that do not contain any terms proportional to 
v', . . . , v (k ~ 2 \ Therefore, in analogy with the above discus- 
sion of the k = 3 and k = 4 cases, the WENO scheme that 
uses stencils of length k can be shown to be able to handle 
the maxima without reducing the reconstruction order if at 
least one of the 1st, 2nd, . . . , (k — 2)th derivatives at x = Xi 
is nonzero. 

Note that for k ^ 3 the smoothness indicators are 
the same for all reconstruction types: average to interface, 
{v)i — > Wi+1/2, cell centre to interface, Vi — > fi+1/2, average 
to centre, (v)i — > v iy and centre to average, Vi — > (v)i, recon- 
structions and are given by equations (|A18[) . This is because 
for k ^ 3 the difference between cell centre and cell average 
values is the same for all points in the stencil. However, for 
larger k, this difference starts to vary from point to point 
within the stencil, and for k ^ 4 the smoothness indica- 
tors become different for these types of reconstruction. So, 
the 4th order smoothness indicators presented here for cell 
centre to interface, Vi — > v i+1 / 2 , reconstruc tion are differ- 
ent fro m the smoothness indicators given bv lBalsara fc Shul 
(2000) for the average to interface, {v)i — > Vj+i/2, recon- 
struction. 

A7 Algorithm for reducing the stencil size 

WENO schemes only operate on stencils of a fixed length 
and this length is the smallest the stencils can get (e.g., it 
is 3 for the case of WENO-5). There are several cases when 
there may not be large enough smooth region for a stencil of 
a large size to fit in, e.g. (1) for a reconstruction between two 
shocks that are about to interact, (2) in the case of a state 
that has to be 'built-up', and (3) for the gr id cells inside 
unresolved discontinuities l|Harten et al.|[l987h . 

We have developed an algorithm of adaptive reduction 
of the stencil size: the algorithm locates the grid cells inside 
sharp discontin uities in th e flow, and a lower-order WENO-3 
reconstruction ( Shu 1997) is used there. In this algorithm we 
do not rely on the p ressure jumps (as is done i n the PPM 
algorithm, see, e.g., IColella fc Woodward! [1983 ) as indica- 
tors of shocks because in the supersonic regime, where the 
pressure is negligible and, possibly, noisy, this approach may 
lead to excessive reduction of the order of the scheme. 

The information about the discontinuities in the flow is 
actually encoded in the weights (|23p used by the WENO re- 
construction. For the purposes of the stencil reduction algo- 
rithm described below, we compute the unoptimized weights 
that are obtained by setting d r = 1 in (|23p . which, given the 
smoothness indicators, is a fast computation. We want to 
find out if we need to reduce the order of reconstruction 
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for the ith grid cell and we are using the WENO-5 scheme 
with k — 3. Suppose, the weight u>q is very large compared 
to 0J2 +2 , i-e. the flow around the ith cell has a steeper pro- 
file than that near the (i + 2)nd cell. This is an indication 
that the ith cell is located inside a discontinuit y, so we us e 
the lower-order WENO-3 reconstruction there (|Shulll997f) . 
Namely, for ith grid cell we define the ratio 



7L = T 



V-2 



(A36) 



where we choose function T below depending on the re- 
construction type, and define the fraction of the WENO-3 
reconstruction that we use in that grid cell as 



«3 = max 



i TZ TZuiin 1 

' mm ' t? t? ' 

'Vmax 'vmin 



(A37) 



so the fraction of the WENO-5 reconstruction is (f — 03). 
This means that for the ratios of weights 1Z < 7?. m i n we fully 
use the WENO-5 reconstruction procedure, for 1Z > 7?. ma x 
we fully use the WENO-3 one, and linearly transition be- 
tween them for intermediate 1Z. For the cell centre to inter- 
face reconstruction we use J-(x,y) = ma.x(x,y), 7?. m in = 1-3, 
and 7£ ma x = 1.6. This provides reduction in the order of 
reconstruction in unresolved discontinuities whose width is 
smaller than about the stencil size (5 points for WENO- 
5) and avoids reduction in the kinks of the flow (disconti- 
nuities in first derivative). For centre to average and aver- 
age to centre reconstructions we use J-(x,y) = mm(x,y), 
72. m i n = 10, and 7?. max = 15. This allows reduction in unre- 
solved discontinuities and in the kinks of the flow. Lastly, 
the stiff relativistic regime where the effective Lorentz fac- 
tor 7 c ff = ~/(p + u g +p g )/p ^ 10 stresses higher-order meth- 
ods since, unlike with TVD methods, the interpolated value 
can lie outside of the range given by the surrounding val- 
ues. To avoid mild oscillations in the ultrarelativistic test 
in section 15.191 we use the latter reduction method for all 
reconstructions if 7 e g ^ 10. This helps by keeping the inter- 
polations consistent with each other and only makes a mild 
difference near the shock in this single ultrarelativistic test. 

In smooth flows where the first and second derivatives 
do not vary much from grid cell to grid cell WENO-5 main- 
tains high order accuracy: all weights are on the order of 1/3, 
and the reduction does not operate. However, if there is a 
discontinuity in the flow, the weights for the points that are 
near it will vanish for the stencils crossing the discontinuity, 
so a lower order reconstruction will be used for the points 
inside of the discontinuity. Simil ar to what is done in the flat- 
tening algorithm use d by PPM (|Colella fe Woodward1ll984l ; 
iMignone et"alll2005h . we additionally reduce the order near 
shocks but avoid such a reduction near contact discontinu- 
ities. For each grid cell we modify Q3 as follows: 



Q3 = max(a 3 , <Sja3,i-i, SiOi3,i+i) 



(A38) 



where a3,i±i are the unmodified values of as at the grid cells 
adjacent to the ith one and Si is an indicator of the shock 
strength that vanishes for no shock and becomes unity for a 
strong shock: 



Si — max |o, min(l, 45 4 — l)j , 

PiA |u*(l + u t )\ + A \(p g + u s )ii*Mi 



Si = 



[pU*(l + Ut) + {p g + Ug)li*Ut +Pg]i 



(A39) 
(A40) 



Si has the meaning of the relative jump in the energy (— T'f — 
pu l ), eq. 0, and Ag = q l+1 - qi-i. 

We introduce reduced WENO-5 weights, iu r = (1 — 
ctz)ui r , and rewrite the reconstruction result, 

/weno = h + cb/3, (A41) 

where fs is the reconstruction due to WENO-5 with the 
reduced weights, is the reconstructed value due to the 
WENO-3, and &3 = 1 — ~}2 r ui r . We will modify the reduced 
WENO-5 weights in the following paragraph. 

We perform the integration of source terms in the usual 
component-by-component way. However, we have to be more 
careful for average to centre and centre to average recon- 
structions since they operate on the conserved quantities. 
The difference between cell averaged and cell centred values 
cannot be treated in a pure component-by-component way 
because of the nonlinear coupling between the conserved 
quantities. For example, we would like to avoid using dif- 
ferent stencils in average to centre reconstructions for to- 
tal conserved energy and components of the total conserved 
momentum. For the conserved quantities we separately per- 
form the reconstruction of the energy component in the way 
described in the previous paragraph. For each of the other 
components we modify the reduced WENO-5 weights to be 
equal to the minimum of the reduced WENO-5 weights for 
the energy component and ten times the WENO-5 weights 
for that component. This way in smooth regions for the re- 
construction of all conserved quantities we use a common set 
of weights computed for the total conserved energy. Simi- 
larly, for the centre to average reconstruction of fluxes in 
ith direction, we use the same procedure as described above 
for the conserved quantities, but instead of using a com- 
mon set of weights computed for the energy component, we 
use a common set of weights computed for the flux of ith 
component of momentum in the i the direction. 

We make the reconstruction even more robust by re- 
ducing its order in non-smooth features of the flow that 
we refer to as cusps (defined below). Further, if there is 
a cusp in the 7-factor of the flow in a grid cell and cell to 
cell change of 7 is larger than 1%, we do not perform the 
centre to average reconstruction of the fluxes whose values 
depend on the reconstruction of 7 that has the cusp (since 
in this case the fluxes, which are steep functions of 7, are 
not smooth). We now explain what we mean by a cusp. De- 
fine fl +1/2 = f i+ i - f iy f" = f' i+lh - /■_!/,. First, consider 



the case of an increasing function, f[_\u > V^WfW', we later 
will generalize the procedure. We then define a point i to 
be in a cusp if it is located between an inflection point (i.e., 
fi-i > \/e||/|| and /" < — \/s||/||) and a local maximum 
with one-sided derivatives different by no more than 50% by 
absolute value (i.e. either 4\f' i+ i /2 +./t V ,| < |/- + i/ 2 | + |/-_i/ 2 1 

or 4\f' i+3/2 + fl_x h \ < I/-+3/J + \fU/ 2 \)- In S eneraL we de- 
fine the function to contain the cusp at the point i if either 
the above is true or if it becomes true after the flip of the 
function sign, grid direction, or both. For consistency, we 
choose e to be the same as the one in the WENO weights 
computation procedure (|A15|) . (|A16I) . (J23J), QZty . 



A8 Failsafe integration 

As in oth er codes (|Zhang fc MacFadvenl 120061 : 

iGammie et ail 120031 ). we implement certain safety fea- 
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tures that enable our code to succeed in especially stiff 
regimes. In relativistic flows, there are states of conserved 
quantities for which there exists no corresponding set 
of primitive quantities (|10[) at all. Further, a change in 
conserved quantities of less than 1% may correspond to 
an infinite change in primitive quantities. Therefore, any 
reconstruction operations directly performed on the con- 
served quantities or on the fluxes may in principle lead to 
unphysical states or instabilities. This is the reason why we 
carry out the crucial cell centre to interface reconstruction 
|Stcp i| (see section 13. 6[) on the primitive quantities instead 
of conserved quantities or fluxes: this guarantees that we 
always obtain a physically meaningful state at the interface 
and, therefore, a physically meaningful flux. However, the 
update to the conserved quantities due to this flux may 
sometimes lead to an unphysical state with a negative 
internal energy or density, or, as explained above, the 
corresponding state of primitive quantities may not at all 
exist. Here we describe the procedure that we follow in such 
cases in order to make the integration failsafe. 

In our scheme the de-averaging and averaging recon- 
structions, [StepTn] and [Step^Ii] (see section [3~6ll . change the 
representation of conserved quantities but do not alter the 
locations of those. We avoid using the average to centre re- 
construction if it makes a large difference for the total energy 
density in the fluid frame, (p + u g ) in the hydrodynamic case 
we consider in this paper. In particular, we make sure that 
the primitive quantities p c , that correspond to cell centred 
conserved quantities, are not very different from p a , that 
correspond to cell averaged conserved quantities: we use p c 
if its difference from p a in the value of energy density in the 
fluid frame is smaller than 5%, use p a if that difference is 
larger than 10%, and linearly transition between these two 
values for an intermediate difference. If no primitive state 
p c has been found by the inversion, we use p a . This locally 
lowers the order of the scheme in stiff regimes where opera- 
tions on the conserved quantities are not safe, so there our 
scheme becomes equivalent to the schemes that do not dis- 
tinguish between the average and point values. Note, that 
this procedure does not affect the asymptotic order of con- 
vergence of the scheme since the difference between the cell 
averaged and cell centred values is 0(h 2 ). 

Occasionally, during the evolution negative rest-mass 
density or internal energy might occur. If a negative rest- 
mass density (internal energy) implies an imaginary sound 
speed, we set the sound speed to unity (zero) and continue 
the evolution. This adjustment of the sound speed only af- 
fects the diffusive flux used in the approximate HLL Rie- 
mann solver. Most of such negative rest-mass density and 
internal energy occurrences are due to inaccuracy of trial 
Runge-Kutta time steps and do not appear on the final 
Runge-Kutta time steps. Such behaviour is acceptable since 
at the end we get a high order accurate answer, and it does 
not matter that some of the intermediate steps that lead to 
this answer were unphysical. For the same reason we allow 
negative internal energy /baryon densities to occur at the 
final Runge-Kutta time steps. The idea is the same: to tem- 
porarily allow the solution to go out of the physical states 
space with the hope of it finally reaching a well-defined ac- 
curate state. For instance, negative internal energy/baryon 
densities may occur in front of a strong shock, however, we 
find that once the shock passes through, all of the states that 



were temporarily unphysical are well-defined. Lastly, we im- 
plement a minor diffusive correction to the internal energy 
under special circumstances in order the improve the fidelity 
of the internal energy near shocks (as for a few cells near the 
shock in the caustic test in section 15. 3[) . If the internal en- 
ergy is negative but all surrounding values are larger, then 
the internal energy is considered to be an isolated failure 
and its value is chosen to be the smallest of the surrounding 
values in a way that leads to a symmetric result in multiple 
dimensions. 

Finally, we apply some diffusive corrections to the prim- 
itive quantities used to obtain the flux. In rare cases that no 
state of primitive quantities can be found that corresponds 
to the state of the conserved quantities (i.e. an inversion 
from conserved to primitive quantities does not exist) we use 
the average of the existing surrounding states of the prim- 
itive quantities. This situation occurred only twice in the 
two-dimensional relativistic Riemann problem (section l5.2UI) 
in the highly relativistic part of the flow (7 ~ 25) and did 
not occur in any other test problems we have performed. 
If for all surrounding grid cells no primitive states can be 
found, we use the average of states of primitive quantities in 
the surrounding grid cells from the previous time step. 

Note that the above features only alter the state of 
primitive quantities that are used in order to compute in- 
terface fluxes, therefore, the scheme remains conservative 
since the values of the conserved quantities are affected only 
through the fluxes. 

A9 Catastrophic cancellation in nonrelativistic 
problems 

In order to be able to successfully study both highly 
relativistic and nonrelativistic problems and avoid being 
severely limited by finite numerical precisiorjf], we have to 
use the a special technique for performing the rest-mass 
subtraction in © and l[8jl. For instance, in the calculation 
of conserved quantities from the primitive ones, given a 4- 
velocity tt M , we need to find the kinetic energy, i.e. the differ- 
ence between the total bulk motion energy and the rest-mass 
energy in the coordinate basis frame, 

pu(-u t - 1). (A42) 

A similar procedure is followed in lAlov et al.l (|l999T ) for the 
special relativistic equations of motion. Further, the inverse 
process of computing the primitive quantities from con- 
served quantities involves computing a similar expression 

p 7 (7-l). (A43) 

For nonrelativistic flows the value of 7 may get so close 
to unity that within machine precision its numerical value 
is 1.0. This is why we cannot simply plug in such a value 
for 7 in the above formula to compute the kinetic energy. 
For double precision, this happens for 3-velocity values of 
v < 10 -8 . We alleviate the above problem by rewriting part 

9 Even though we normally use the code at double precision, we 
have made it capable of performing computations at long double 
precision. For this, we use the long double versions of exponential, 
trigonometric, etc., functions from Cephes Math Library Release 
2.7 and use the -long_double option for the Intel compiler. 
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of expression (|A43[) for kinetic energy with the help of the 
velocity information: 

1 - yi-v 2 



= - y - » = J" = j__ / A44 n 

v / r^ 2 i + ^r~^2 i +7 ' ^ i 

where the final expression does not contains catastrophic 
cancellations and due to the particular choice of frame 
the Lor entz factor has th e conventional form 7 = (1 — 
v 2 )' 1 ' 2 dNoble et alj|2006h . 

In equation (|A42|) we can similarly rewrite the problem- 
atic part as follows: 



2 2 
7 v 



Ut 



1 = — <7«u' + 20- 



- gtt i 2 v 2 



, , (A45) 
1+7 1+7 

where 7 = —guU* and £> 2 = 2guv l — 2(f> + gijV^v^ play the 
roles of the 7-factor and the square of the velocity, <f> = 
— (1 + gtt)/ 2 is the gravitational potential and is small com- 
pared to unity for problems involving nonrelativistic gravity, 
so for such problems we separately store its value. Simi- 
larly, we separately compute and store the value of (g — 1), 
which appears when converting the total energy with rest- 
mass s ubtracted from th e lab frame to the normal ob server 
frame i|Noble et al.ll2006l ; iMignone fc McKinney||200^ . 

The general expression ()A45|) avoids catastrophic can- 
cellations both for relativistic problems and for problems in- 
volving nonrelativistic velocities and gravitational fields. We 
use this expression when computing the conserved quanti- 
ties from the primitive ones as well as for the inverse process. 
Since expression (|A45[) is written in an arbitrary frame, it in- 
volves two additional terms as compared to (|A44|) : the first 
one appears for a metric that has space-time mixing, the 
second one is due to a nonzero gravitational potential, and 
the third one corresponds to (|A44|) . 

A10 Catastrophic cancellation in ultrarelativistic 
problems 

For the ultrarelativis tic regime one can show that the inver- 
sion method we use (IMignone fc McKinnevll2007l ) that con- 
verts conserved to primitive quantities leads to well-defined 
relative errors in the primitive quantities for a machine ac- 
curate set of conserved quantities. In the ultrarelativistic hy- 
drodynamic case, for a conserved lab-frame energy density 
(E), conserved mass density (D), and conserved momentum 
density (P a ), the Lorentz factor is given by 

771 

j= — =, (A46) 

which can be used to obtain the rest-mass density 

D 



P = 



7 



the relative 4-velocity 



~ a ~a 

u = 71; =7 



pa 

~E 



and the quantity \ — u g + Pg , 



_ E_ _ D 

7 7 



(A47) 



(A48) 



(A49) 



which with an equation of state can be used to determine 
both u g and p g . For this inversion we have neglected the 
term in the energy equation proportional to p g compared to 



7 2 (p + X) an d below we state how this limits the remaining 
arguments. 

Given E, P, and D with an accuracy limited by a ma- 
chine error of e ma chm C , one can show that the relative errors 
in primitive quantities are given by 



d'y 
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(A50) 
(A51) 
(A52) 
(A53) 



These error estimates are valid as long as 

— » (^r-) 7 • ( A54 ) 

This implies a smaller than order unity error for the primi- 
tive quantities when 

= ,-1/2 



p(7/7max) 2 
1 - (7/7max) 2 ' 



(A55) 
(A56) 



and the error in p only sets the maximum allowed 7. For 
example, using double precision on a 32-bit machine gives 
conserved quantities that have at best a relative error of 
^machine ~ 2.2 x 10" 16 giving 



7max = 6.7 X 10 7 . 



(A57) 



For 



sampl e, 



th e ultrarelativistic test discussed 
in lAlov et all l|l999l ) (see section EH91 with 7 = 2.24 x 10 5 
and p = 7.63 x 10 -6 for F = 4/3 gives x ^ Xmin, so their 
test is at the limit of resolving the pre-shock pressure for 
double precision. The post-shock region will be resolved as 
long as the pre-shock Lorentz factor is resolved. 

Finally, an additional source of error can be incurred 
when an iterative inversion method uses the expression 

,.2\-l/2 



7 = (1- v y 



as compared to 



(A58) 



(A59) 



since repeated use of the prior expression leads to a cumu- 
lated catastrophic errors not present when using the latter 
expression. For this reason, in p ractice the GRMHD itera- 
tive solver of iNoble et~all l|2006l) is limited to 7 < 10 5 when 
using double precision. 
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